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There are some things which cannot 
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Preface 



The vast majority of this book comes from lecture notes I have been typing 
up over the years for a could on differential equations with boundary value 
problems at the USNA. Though the USNA is a government institution and 
official work-related writing is in the public domain, I typed and polished 
so much of this at home during the night and weekends that 1 feel 1 have 
the right to claim copyright over this work. The DE course has used various 
editions of the following three books (in order of most common use to least 
common use) at various times: 

• Dennis G. Zill and Michael R. Cullen, Differential equations with 
Boundary Value Problems, 6th ed., Brooks/Cole, 2005. 

• R. Nagle, E. Saff, and A. Snider, Fundamentals of Differential 
Equations and Boundary Value Problems, 4th ed., Addison/Wesley, 
2003. 

• W. Boyce and R. DiPrima, Elementary Differential Equations and 
Boundary Value Problems, 8th edition, John Wiley and Sons, 2005. 

You may see some similarities but, for the most part, I have taught things a 
bit differently and tried to impart this in these notes. Time will tell if there 
are any improvements. 

A new feature to this book is the fact that every section has at least one 
SAGE exercise. SAGE is FOSS (free and open source software), available on 
the most common computer platforms. Royalties for the sales of this book 
(if it ever makes it's way to a publisher) will go to further development of 
SAGE. 

This book is free and open source. It is licensed under the Attribution- 
ShareAlike Creative Commons license, http: // creativecommons . org/ licens 
by-sa/3. 0/ , or the Gnu Free Documentation License (GFDL), http: 
//www. gnu. org/ copyleft/ fdl . html, at your choice. 
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Intro... 

If people do not believe that mathematics is simple, it is only 
because they do not reahze how comphcated hfe is. 
- John von Neumann 



To be written 



Chapter 1 

First order differential 
equations 

1.1 Introduction to DEs 

But there is another reason for the high repute of mathe- 
matics: it is mathematics that offers the exact natural sciences 
a certain measure of security which, without mathematics, they 
could not attain. 

- Albert Einstein 



Motivation 

Roughly speaking, a differential equation is an equation involving the deriva- 
tives of one or more unknown functions. 

In calculus (differential, integral and vector), you've studied ways of ana- 
lyzing functions. You might even have been convinced that functions you 
meet in applications arise naturally from physical principles. As we shall 
see, differential equations arise naturally from general physical principles. In 
many cases, the functions you met in calculus in applications to physics were 
actually solutions to a "natural" differential equation. 

Example 1.1.1. Consider a falling body of mass m on which exactly 3 forces 
act: 



1 



• gravitation, Fg^av, 

• air resistance, Fres, 

• an external force, Fext — fif), where fit) is some given function. 











mass m 




F 

^ grav 



Let x{t) denote the distance fallen from some fixed initial position. The 
velocity is denoted by v — x' and the acceleration by a — x" . We choose 
an orientation so that downwards is positive. In this case, Fg^av — mg, 
where g > is the gravitational constant. We assume that air resistance is 
proportional to velocity (a common assumption in physics), and write Fres = 
—kv = —kx', where k > is a 'friction constant". The total force, Ftotai, is 
by hypothesis, 

F'total — F'grav ~l~ F'res ~l~ F^xt: 

and, by Newton's 2nd Law^, 

Ftotai ^ma^ mx". 
Putting these together, we have 

mx" — ma — mg — kx' + f{t), 

or 

mx" + mx' = fit) + mg. 

This is a differential equation in x = x{t). It may also be rewritten as a 
differential equation in v — v{t) — x'{t) as 

^ "Force equals mass times acceleration." http://en.wikipedia.org/wiki/Newtons_ 
law 



mv' + kv = fit) + mg. 

This is an example of a "first order differential equation in v ", which means 
that at most first order derivatives of the unknown function v = v{t) occur. 

In fact, you have probably seen solutions to this in your calculus classes, at 
least when f{t) — and k — 0. In that case, v'{t) — g and so v{t) — J gdt — 
gt + C. Here the constant of integration C represents the initial velocity. 

Differential equations occur in other areas as well: weather prediction (more 
generally, fluid-flow dynamics), electrical circuits, the heat of a homogeneous 
wire, and many others (see the table below). They even arise in problems 
on Wall Street: the Black-Scholes equation is a PDE which models the pric- 
ing of derivatives [BS-intro]. Learning to solve differential equations helps 
understand the behaviour of phenomenon present in these problems. 



phenomenon 


description of DE 


weather 


Navicr-Stokes equation [NS-intro] 
a non-linear vector-valued higher-order PDE 


falling body 


1st order linear ODE 


motion of a mass attached 
to a spring 


Hookc's spring equation 
2nd order linear ODE [H-intro] 


motion of a plucked guitar string 


Wave equation 
2nd order linear PDE [W-intro] 


Battle of Trafalger 


Lanchester's equations 
system of 2 1st order DEs [L-intro], [M-intro], [N-intro] 


cooling cup of coffee 
in a room 


Newton's Law of Cooling 
1st order linear ODE 


population growth 


logistic equation 
non-linear, separable, 1st order ODE 



Undefined terms and notation will be defined below, except for the equations 
themselves. For those, see the references or wait until later sections when 
they will be introduced^. 



Basic Concepts: 

Here are some of the concepts to be introduced below: 

^Except for the Navier-Stokes equation, which is more complicated and might take us 
too far afield. 



dependent variable(s), 



• independent variable(s), 

• ODEs, 

• PDEs, 

• order, 

• linearity, 

• solution. 

It is really best to learn these concepts using examples. However, here are 
the general definitions anyway, with examples to follow. 

The term "differential equation" is sometimes abbreviated DE, for brevity. 

Dependent /independent vairiables: Put simply, a differential equation 
is an equation involving derivatives of one of more unknown functions. The 
variables you arc differentiating with respect to are the independent vari- 
ables of the DE. The variables (the "unknown functions") you are differenti- 
ating are the dependent VEiriables of the DE. Other variables which might 
occur in the DE are sometimes called "parameters" . 

ODE/PDE: If none of the derivatives which occur in the DE are partial 
derivatives (for example, if the dependent variable/unknown function is a 
function of a single variable) then the DE is called an ordinary differential 
equation of PDE. If some of the derivatives which occur in the DE are 
partial derivatives then the DE is a peirtial differential equation or PDE. 

Order: The highest total number of derivatives you have to take in the 
DE is it's order. 

Linearity: This can be described in a few different ways. First of all, a DE 
is linear if the only operations you perform on its terms are combinations of 
the following: 

• differentiation with respect to independent variable(s), 

• multiplication by a function of the independent variable(s). 



Another way to define linearity is as follows. A linear ODE having inde- 
pendent variable t and the dependent variable is y is an ODE of the form 



ao(t)|/(") + ... + a„_i(t)i/' + anit)y = f{t), 
for some given functions ao(t), . . . , an(^), and fit). Here 

denotes the n-th derivative of y = y{t) with respect to t. The terms ao(t), 
. . . , an{t) are called the coefficients of the DE and we will call the term 
f{t) the non- homogeneous term or the forcing function. (In physical 
applications, this term usually represents an external force acting on the 
system. For instance, in the example above it represents the gravitational 
force, mg.) 

Solution: An explicit solution to a DE having independent variable t and 
the dependent variable is x is simple a function x{t) for which the DE is true 
for all values of t. 

Here are some examples: 



Example 1.1.2. Here is a table of examples. As an exercise, determine 
which of the following are ODEs and which are PDEs. 



DE 


indep vars 


dep vars 


order 


linear? 


mx" + kx' = mg 
falling body 


t 


X 


2 


yes 


mv' + kv = mg 
falling body 


t 


V 


1 


yes 


l.d'-'u Ou 
'^dx'^ ^ dl. 

heat equation 


t, x 


u 


2 


yes 


mx" + bx' + kx = f(t) 
spring equation 


t 


X 


2 


yes 


P' = fc(l-^)P 
logistic population equation 


t 


p 


1 


no 


wave equation 


t, X 


u 


2 


yes 


k(T ^room) 

Newton's Law of Cooling 


t 


T 


1 


yes 


x' = -Ay, y' = -Bx, 
Lanchester's equations 


t 


X, y 


1 


yes 



Remark; Note that in many of these examples, the symbol used for the 
independent variable is not made explicit. For example, we are writing x' 
when we really mean x'{t) — This is very common shorthand notation 
and, in this situation, we shall usually use t as the independent variable 
whenever possible. 

Example 1.1.3. Recall a linear ODE having independent variable t and the 
dependent variable is y is an ODE of the form 

ao{t)y^''^ + ... + an-i(t)y' + an(t)y = f(t), 

for some given functions ao{t), . . . , an{t), and f{t). The order of this DE is 
n. In particular, a linear 1st order ODE having independent variable t and 
the dependent variable is y is an ODE of the form 

ao{t)y' + ai{t)y = f{t), 

for some aQ{t), ai{t), and f{t). We can divide both sides of this equation by 
the leading coefficient aoit) without changing the solution y to this DE. Let's 
do that and rename the terms: 

y' +p{t)y^q{t), 

where p{t) — ai{t)/ao{t) and q{t) — f{t)/ao{t). Every linear 1st order ODE 
can be put into this form, for some p and q. For example, the falling body 
equation mv'+kv — f{t)+mg has this form after dividing by m and renaming 
V as y. 

What does a differential equation like mx" + kx' — mg or P' = k{l — 
or = ^ really mean? In mx" + kx' = mg, m and k and g are given 
constants. The only things that can vary are t and the unknown function 
X = x(t). 

Example 1.1.4. To be specific, let's consider x' + x = 1. This means for all 
t, x'(t) + x(t) = 1. In other words, a solution x{t) is a function which, when 
added to its derivative you always get the constant 1. How many functions 
are there with that property? Try guessing a few "random" functions: 

• Guess x{t) = sin(t). Compute (sin(t))' + sin(i) = cos(i) + sin(i) = 
V2sin(t + f ). x'{t) + x{t) = 1 is false. 



• Guess x{t) — exp(i) = e*. Compute (e*)' + e* = 2e*. x'{t) + x{t) — 1 is 
false. 

• Guess x{t) = exp(i) = Compute {t'^y + t'^ ^2t + t'^. x'{t)+x{t) = 1 
is false. 

• Guess x{t) — exp(— t) = e~*. Compute (e~*)' + e~* = 0. + = 1 
is false. 

• Guess x{t) = exp(t) = 1. Compute (1)' + 1 = + 1 = 1. x'{t) + x{t) = 1 
is true. 

We finally found a solution by considering the constant function x{t) = 1. 
Here a way of doing this kind of computation with the aid of the computer 
algebra system SAGE ; 

SAGE 

sage : t = var { ' t' ) 

sage: de = lambda x: diff (x,t) + x - 1 

sage : de (sin (t) ) 

sin(t) + cos(t) - 1 

sage : de (exp (t) ) 

2*e"t - 1 

sage : de (t " 2 ) 

f2 + 2*t - 1 

sage: de(exp(-t)) 

-1 

sage : de ( 1 ) 




Note we have rewritten x' + x — 1 as x' + x — 1 — and then plugged various 
functions for x to see if we get or not. 

Obviously, we want a more systematic method for solving such equations 
than guessing all the types of functions we know one-by-one. We will get to 
those methods in time. First, we need some more terminology. 

IVP: A first order initial value problem (abbreviated IVP) is a problem 
of the form 



x' = fit, x), x{a) = c, 

where f(t, x) is a given function of two variables, and a, c are given constants. 
The equation x{a) = c is the initial condition. 

Under mild conditions of /, an IVP has a solution x — x{t) which is unique. 
This means that if / and a are fixed but c is a parameter then the solution 
X — x{t) will depend on c. This is stated more precisely in the following 
result. 

Theorem 1.1.1. (Existence and uniqueness) Fix a point (to, xq) in the plane. 
Let f{t, x) be a function oft and x for which both f{t, x) and fx{t, x) — ^^g^^^ 
are continuous on some rectangle 

a < t < b, c < X < d, 

in the plane. Here a, b, c, d are any numbers for which a < to < b and 
c < xq < d. Then there is an h > and a unique solution x — x{t) for which 

x' — f{t, x), for all t e (to — h, tQ + h), 

and x{to) = Xq. 

This is proven in §2.8 of Boyce and DiPrima [BD-intro], but we shall not 
prove this here. In most cases we shall run across, it is easier to construct 
the solution than to prove this general theorem. 

Example 1.1.5. Let us try to solve 

x' + x^l, x(0) = 1. 

The solutions to the DE x' -\- x — 1 which we "guessed at" in the previous 
example, x{t) = 1, satisfies this IVP. 

Here a way of finding this slution with the aid of the computer algebra 
system SAGE ; 

SAGE 

sage : t = var { ' t' ) 

sage: x = f unction (' x' , t) 

sage: de = lambda y: diff(y,t) + y - 1 

sage: desolve_laplace (de (x (t) ) , ["t", "x"] ,[0,1]) 



(The command desolve_laplace is a DE solver in SAGE which uses a special 
method involving Laplace transforms which we will learn later.) Just as an 
illustration, let's try another example. Let us try to solve 

x' + x = l, x(0) = 2. 
The SAGE commands are similar: 

SAGE 

sage : t = var { ' t' ) 

sage: x = f unction {' x' , t) 

sage: de = lambda y: diff(y,t) + y - 1 

sage: desolve_laplace (de (x (t) ) , [ "t", "x"] , [0,2]) 

' %e'-t+l' 

age: solnx = lambda s: RR (eval { soln . replace {""","**") . 

replace {"%", "") . replace ( "t " , str (s) ) ) ) 

sage: solnx{3) 

1 . 04978706836786 

sage: P = plot ( solnx, , 5 ) 

sage: show{P) 




Figure 1.1: Solution to IVP x' + x = l, x(0) = 2. 



Exercise: Verify the, for any constant c, the function x{t) = 1 + ce~* solves 
x' + X = 1. Find the c for which this function solves the I VP x' + x — 1, 
x{0) = 3.. Solve this (a) by hand, (b) using SAGE. 

1.2 Initial value problems 

A 1-st order initial value problem, or IVP, is simply a 1-st order ODE 
and an initial condition. For example, 

x'{t) + p{t)x{t) = q{t), x{0) = Xq, 

where p{t), q{t) and xq are given. The analog of this for 2nd order linear 
DEs is this: 

a{t)x"{t) + h(t)x'{t) + c(t)x{t) = fit), x{Q) = xo, x\0) = vq, 

where a{t), b{t), c{t), xq, and vq are given. This 2-nd order linear DE and 
initial conditions is an example of a 2-nd order IVP. In general, in an IVP, 
the number of initial conditions must match the order of the DE. 

Example 1.2.1. Consider the 2-nd order DE 

x" + x^ 0. 

(We shall run across this DE many times later. As we will see, it represents 
the displacement of an undamped spring with a unit mass attached. The term 
harmonic oscillator is attached to this situation [0-ivpJ.) Suppose we know 
that the general solution to this DE is 

x{t) — Ci cos(i) -I- 02 sin(t), 

for any constants ci, c^,. This means every solution to the DE must he of this 
form. (If you don't believe this, you can at least check it it is a solution by 
computing x" it) + x{t) and verifying that the terms cancel, as in the following 
SAGE example. Later, we see how to derive this solution.) Note that there 
are two degrees of freedom (the constants Ci and C2), matching the order of 
the DE. 

SAGE 

sage : t = var ( ' t' ) 



sage: cl = var{'cl') 

sage: c2 = var{'c2') 

sage: de = lambda x: diff(x,t,t) + x 

sage: de(cl*cos{t) + c2*sin(t)) 



sage: x = f unction {' x' , t) 

sage: soln = desolve_laplace (de (x (t ) ) , [ "t " , "x" ] , [ , , 1 ] ) 

sage: soln 
' sin (t) ' 

sage: solnx = lambda s: RR (eval { soln . replace { "t "," s ")) ) 

sage: P = plot { solnx, , 2 *pi ) 

sage: show{P) 



This is displayed below. 
Now, to solve the IVP 

x" + x = 0, x{0) = 0, a;'(0) = 1. 

the problem is to solve for C\ and for which the x{t) satisfies the initial 
conditions. The two degrees of freedom in the general solution matching the 
number of initial conditions in the IVP. Plugging t = into x{t) and x'{t), 
we obtain 

= x(0) = Cl cos(O) + C2 sin(O) = ci, 1 = x'(0) = — Ci sin(O) + C2 cos(O) = C2. 
Therefore, ci = 0, C2 = 1 and x{t) = sin(t) is the unique solution to the IVP. 




Figure 1.2: Solution to IVP x" + x = 0, x(0) = 0, x'(0) = 1. 
Here you see the solution oscillates, as t gets larger. 



Another example, 
Example 1.2.2. Consider the 2-nd order DE 

x" + 4x' + 4a; = 0. 

(We shall run across this DE many times later as well. As we will see, it 
represents the displacement of a critially damped spring with a unit mass 
attached.) Suppose we know that the general solution to this DE is 

x{t) = ciexp{-2t) + C2texp{-2t) = Cie~^* + C2te~^\ 

for any constants c\, C2. This means every solution to the DE must he of 

this form. (Again, you can at least check it is a solution by computing x"{t), 
4a;' (t), 4a; (t), adding them up and verifying that the terms cancel, as in the 
following SAGE example.) 

SAGE 

sage : t = var ( ' t' ) 
sage: cl = var('cl') 
sage: c2 = var('c2') 

sage: de = lambda x: diff{x,t,t) + 4*diff(x,t) + 4*x 

sage: de {cl*exp {-2*t) + c2*t*exp (-2*t) ) 

4* {c2*t*e" {-2*t) + cl*e" (-2*t) ) + 4* {-2*c2*t*e" {-2*t) 

+ c2*e'(-2*t) - 2*cl*e' (-2*t) ) + 4*c2*t*e' (-2*t) 

- 4*c2*e' {-2*t) + 4*cl*e" {-2*t) 

sage: de (cl*exp (-2*t) + c2*t*exp (-2*t) ). expand ( ) 



sage: desolve_laplace (de (x (t) ) , ["t", "x"] , [0,0,1]) 
't*%e''- (2*t) ' 

sage: P = plot (t*exp (-2*t) , 0, pi) 
sage: show{P) 



The plot is displayed below. 
Now, to solve the IVP 

x" + 4a;' + 4a; = 0, a;(0) = 0, a;'(0) = 1. 

we solve for c\ and C2 using the initial conditions. Plugging t 
and x'{t), we obtain 



— into x{t) 



= a;(0) = ci exp(O) + Ca ■ ■ exp(O) = Ci, 

1 = x'(0) = ci exp(O) + C2 exp(O) — 2c2 ■ ■ exp(O) = ci + C2. 

Therefore, ci = 0, ci + C2 = 1 one? so a;(t) = t exp(— 2t) zs the unique solution 
to the IVP. Here you see the solution tends to 0, as t gets larger. 



Figure 1.3: Solution to IVP x" + Ax' + 4x = 0, x(0) = 0, x'(0) = 1. 

Suppose, for the sake of argument, that I lied to you and told you the 
general solution to this DE is 

x(t) = ciexp{-2t) + C2exp{-2t) = ci(e~^* + C2e~^*), 

for any constants c\, C2- (In other words, the "extra t factor" is missing.) 
Now, if you try to solve for the constant c\ and C2 using the initial conditions 
x(0) = 0, x'(0) = 1 you will get the equations 



These equations are impossible to solve! You see from this that you must 
have a correct general solution to insure that you can solve your IVP. 

One more quick example. 

Example 1.2.3. Consider the 2-nd order DE 




-1 



-0.5 



-0.75 



Cl + C2 = 
-2Ci - 2C2 = 1. 



x" -x = 0. 

Suppose we know that the general solution to this DE is 

x{t) = Ciexpit) + C2exp{—t) = Cie^* + C2e~*, 

for any constants Ci, 02- (Again, you can check it is a solution.) 
The solution to the IVP 

x" -x = 0, x{Q) = 0, x\Q) = 1, 

is x(t) = (^'^'^ can solve for Ci and C2 yourself, as in the examples 

above.) This particular function is also called a hyperbolic cosine func- 
tion, denoted cosh.{t) . 

The hyperbolic trig functions have many properties analogous to the usual 
trig functions and arise in many areas of applications [H-ivpJ. For example, 
cosh(i) represents a catenary or hanging cable [C-ivp]. 

SAGE 

sage : t = var ( ' t' ) 

sage: cl = var('cl') 

sage: c2 = var{'c2') 

sage: de = lambda x: diff {x,t,t) - x 

sage: de (cl*exp (-t) + c2*exp(-t)) 


sage: desolve_laplace (de (x (t) ) , ["t", "x"] , [0,0,1]) 

' %e"t/2-%e"-t/2' 

sage: P = plot (e't/2-e' (-t) /2, 0, 3) 

sage: show(P) 



Here you see the solution tends to infinity, as t gets larger. 
Exercise: The general solution to the faUing body problem 

mv' + kv — mg, 

is v(t) = ^ + ce^^*/™". If f (0) = vq, solve for c in terms of Vq. Take 
m = k = Vo = l, g = 9.8 and use SAGE to plot v{t) for < t < 1. 



Figure 1.4: Solution to IVP x" - x = 0, x(0) = 0, x'(0) = 1. 



1.3 First order ODEs - separable and linear 
cases 

Separable DEs: 
We know how to solve any ODE of the form 

y' = f{t). 

at least in principle - just integrate both sides^. For a more 
general type of ODE, such as 

y = f{t,y). 

this fails. For instance, if y' = t + y then integrating both sides 
gives y{t) = J^dt = J y' dt = Jt + ydt = Jtdt + J y{t) dt = 

•^Recall y' really denotes so by the fundamental theorem of calculus, y = / ^ = 
J y' dt — J f{t) dt — F{t) + c, where F is the "anti-derivative" of / and c is a constant of 
integration. 



y{t) dt. So, we have only succeeded in writing y{t) in terms 
of its integral. Not helpful. 

However, there is a class of ODEs where this idea works, with 
some slight modification. If the ODE has the form 

' = %)■ 

then it is called separable^. 

To solve a separable ODE: 

(1) write the ODE (1.1) as | = f||, 

(2) "separate" the fs and the y's: 

h{y) dy = g{t) dt, 

(3) integrate both sides: 



Jh{y)dy = Jg{f)df + C 



I've added a "+C" to emphasize that a constant of inte- 
gration must be included in your anwser (but only on one 
side of the equation). 

The answer obtained in this manner is called an "implicit so- 
lution" of (1.1) since it expresses y implicitly as a function of 
t. 

Example 1.3.1. Are the following ODEs separable? If so, solve 
them. 



*It particular, any separable DE must be first order. 



(a) if + y^)y' = -2ty, 
(h) y' = -x/y, yiO) = -1, 

(c) T' = k ■ {T — Troom), where k < and Tj-oom o,re constants, 

(d) ax' -\-hx = c, where a 7^ 0, 6 7^ 0; and c are constants 

(e) ax' -\-hx = c, where a 0, b, are constants and c = c{t) is 
not a constant. 

(f) y' = {y-l){y + l),y{0) = 2. 
(9)y' = y' + i> yio) = i. 

Solutions; 

(a) not separable, 

(b) ydy = —xdx, so = —x'^/2 + c, so x'^ -\- y'^ = 2c. This 
is the general solution (note it does not give y explicitly as 
a function of x, you will have to solve for y algebraically to 
get that). The initial conditions say when x = 0, y = 1, 
so 2c = 0^ + 1^ = \, which gives c = 1/2. Therefore, 

+ y^ = 1, which is a circle. That is not a function 
so cannot be the solution we want. The solution is either 
y = yjl — x"^ or y = — Vl — x'^, but which one? Since 
y{0) = —1 (note the minus sign) it must bey = —yjl — x^. 

(c) 7jT^ — = kdt, so In IT — Troom\ = kt + c (some constant 

^ room 

c), so T — Troom = Cc^* (somc constant C), so T = T{t) = 

room ^ c . 



(d) f = {c-bx)/a = -lix-l), so^ = -Idt, so\n\x-l\ = 

b 

— H-\-C, where C is a constant of integration. This is the 
implicit general solution of the DE. The explicit general 
solution is X = Be~»*, where B is a constant. 

The explicit solution is easy find using SAGE ; 

SAGE 



sage : a = var ( ' a' ) 

sage : b = var { ' b' ) 

sage : c = var { ' c' ) 

sage : t = var ( ' t' ) 

sage: x = f unction (' x' , t) 

sage: de = lambda y: a*diff{y,t) + b*y - c 

sage : desolve_laplace (de (x (t ) ) , [ " t " , "x" ] ) 
' c/b- (a*c-x (0) *a*b) *%e''- (b*t/a) / {a*b) ' 



(e) If c = c{t) is not constant then ax'+bx = c is not separable. 

(f) {y-'l+D = - 1) - Hy +l))=t + C, where C 
is a constant of integration. This is the "general (implicit) 
solution" of the DE. 

Note: the constant functions y{t) = 1 and y{t) = —1 are 
also solutions to this DE. These solutions cannot be ob- 
tained (in an obvious way) from the general solution. 

The integral is easy to do using SAGE ; 
SAGE 

sage : y = var ( ' y' ) 

sage: integral (1/ { (y-l) * (y+l) ) ,y) 
log(y - 1) /2 - (log(y + 1) /2) 



Now, let's try to get SAGE^o solve for y in terms oft in 
i(ln(y-l)-ln(y+l))=t + C; 

SAGE 

sage: C = var ('C ) 

sage: solve{[log{y - l)/2 - {log{y + l)/2) == t+C],y) 
[log(y + 1) == -2*C + log(y - 1) - 2*t] 



This is not working. Let's try inputting the problem in a 
different form: 

SAGE 

sage : C = var { ' C ) 

sage: solve ([ log { (y - l)/{y + 1)) == 2*t+2*C],y) 
[y (-e'(2*C + 2*t) - l)/(e'(2*C + 2*t) - 1)] 



This is what we want. Now let's assume the initial condi- 
tion y{0) = 2 and solve for C and plot the function. 

SAGE 

sage: solny=lambda t : (-e' (2*C+2*t) -1) / (e' (2*C+2*t) -1) 

sage: solve {[ solny { ) == 2],C) 

[C == log {-1/sqrt (3) ) , C == -log(3)/2] 

sage: C = -log{3)/2 

sage: solny (t) 

(-e'(2*t)/3 - 1) / (e' {2*t) /3 - 1) 
sage: P = plot (solny (t) , 0, 1/2) 
sage: show(P) 



This plot is shown below. The solution has a singularity at 
t = ln{3)/2 = 0.5493.... 




Figure 1.5: Plot of y' = {y - l)(y + 1), y{0) = 2, for < t < 1/2.. 

(9) = dt so arctan(?/) = t + C, where C is a constant of 
integration. The initial condition y{0) = 1 says arctan(l) = 
C, so C = J. Therefore y = tan(t + |) is the solution. 

A special subclass of separable ODEs is the class of automo- 
mous ODEs, which have the form 

y' = f{y), 

where / is a given function (i.e., the slope y only depends on 
the value of the dependent variable y). The cases (c), (d), (f), 
and (g) above are examples. 

Linear 1st order DEs: 

The bottom line is that we want to solve any problem of the 
form 



x' +pit)x = q{t), (1.2) 

where p(t) and q{t) are given functions (which, let's assume, 
aren't too horrible). Every first order linear ODE can be writ- 
ten in this form. Examples of DEs which have this form: Falling 
Body problems, Newton's Law of Cooling problems. Mixing 
problems, certain simple Circuit problems, and so on. 
There are two approaches 

• "the formula" , 

• the method of integrating factors. 

Both lead to the exact same solution. 

"The Formula": The general solution to (1.2) is 

^- ^Jp{t)dt ' ^^-^^ 

where C is a constant. The factor e-^^*^*)*^* is called the inte- 
grating factor and is often denoted by /i. This formula was 
apparently first discovered by Johann Bernoulli [F-lst]. 

Example 1.3.2. Solve 

xy' -\- y = e^. 

We rewrite this as y' + = ^. Now compute /j, = = 
gin(a;) _ formula gives 

fx^dx + C fe'^dx + C e^ + C 
y = = = • 



Here is one way to do this using SAGE ; 



SAGE 



sage : t = var ( ' t' ) 

sage: x = f unction (' x' , t) 

sage: de = lambda y: diff{y,t) + (l/t)*y - exp(t)/t 
sage: desolve (de (x (t) ) , [x, t] ) 
' (%e't+%c) /t' 



"Integrating factor method": Let fi = e-^^^*^*^*. Multiply both 
sides of (1.2) by /i: 

jix' -\- pit) fix = fJ^qit). 
The product rule implies that 



(In response to a question you are probably thinking now: No, 
this is not obvious. This is Bernoulli's very clever idea.) Now 
just integrate both sides. By the fundamental theorem of calcu- 
lus, 



(fix)' = fix' + p{t)fix = fiq{t). 




Dividing both side by fi gives (1.3). 



Exercise: (a) Use SAGE's desolve command to solve 

tx' -\-2x = e^/t. 
(b) Use SAGE to plot the solution toy' = y^ -I, y{0) = -2. 



1.4 Isoclines and direction fields 



Recall from vector calculus the notion of a two-dimensional vec- 

— * — * 

tor field: F{x,y) = {g{x,y),h{x,y)). To plot F, you simply 

— * 

draw the vector F{x,y) at each point {x,y). 

The idea of the direction field (or slope field) associated to 
the first order ODE 

y' = f{x,y), y{a) = c, (1.4) 

is similar. At each point (x, y) you plot a vector having slope 
/(x, y). For example, the vector field plot of F(x, y) = (1, /(x, y)) 
or F{x,y) = (1, /(x, y)) / \Jl + f{x,yY (which is a unit vector). 

A related notion are the isoclines of the ODE. An isocline of 
(1.4) is a level curve of the function z = /(x, y): 

{{x,y) I f{x,y) = m}, 

where the given constant m is called the slope of the isocline. In 
terms of the ODE, this curve represents the collection of points 
at which the solution has slope m. In terms of the direction 
field of the ODE, it represents the collection of points where the 
vectors have slope m. 
How to draw the direction field of (I.4) by hand: 

• Draw several isoclines, making sure to include one which 
contains the point (a, c). (You may want to draw these in 
pencil.) 

• On each isocline, draw "hatch marks" or "arrows" along 
the line each having slope m. 



This is a crude direction field plot. The plot of arrows form 
your direction field. The isoclines, having served their useful- 
ness, can safely be ignored or erased. 

Example 1.4.1. The direction field, with three isoclines, for 

y' = ^x + y-5, 7/(0) = 1, 
is given by the following graph: 
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Figure 1.6: Plot of = 5x + ?/ — 5, |/(0) = 1, for — 1 < x < 1. 



The isoclines are the curves (coincidentally, lines) of the form 
^x + y — b = m. (They are green hands in the above plot.) These 
are lines of slope —5, not to he confused with the fact that it 
represents an isocline of slope m. 

The above example can be solved explicitly. (Indeed, y = 
—5x + e^ solves y' = 5x + y — 5, y{0) = 1.) In the next example, 
such an explicit solution is (as far as I know), not possible. 
Therefore, a numerical approximation plays a more important 
role. 



Example 1.4.2. The direction field, with three isoclines, for 

y' = x^ + y\ ?/(0) = 3/2, 
is given by the following graph: 
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Figure 1.7: Direction field and solution plot of y' = + y"^, y{0) = 3/2, for 
-3 < X < 3. 

The isoclines are the concentric circles x'^ + y'^ = m. (They are 
green in the above plot.) 

The plot above was obtained using SAGE 's interface with Max- 
ima, and the plotting package Openmath (SAGE includes both 
Maxima and Openmath). : 

SAGE 

sage : maxima . eval { ' load { "plotdf " ) ' ) 

sage: maxima . eval {' plotdf (x" 2+y" 2 , [tra jectory_at, 0, ] , 

[x,-3,3], [y,-3,3])') 



This gave the above plot. (Note: the plotdf command goes on 
one line; for typographical reasons, it was split in two.) 



There is also a way to draw these direction fields using SAGE . 
SAGE 

sage: pts = [ (-2+1/5, -2+ j/5) for i in range{20) 

for j in range{20)] # square [-2, 2] x [-2, 2] 
sage: f = lambda p : p [ ] " 2+p [ 1 ] " 2 

sage: arrows = [arrow{p, (p [ ] +0 . 02 , p [ 1 ] + { . 02 ) *f (p) ) , 

width=l/100, rgbcolor= ( , , 1 ) ) for p in pts] 
sage: show { sum (arrows ) ) 



This gives the plot below. 
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Figure 1.8: Direction field for y' = x"^ + y'^, y{0) = 3/2, for —2<x<2. 
Exercise: Using SAGE , plot the direction field for y' = — y^. 



1.5 Numerical solutions - Euler's method and 
improved Euler's method 



Read Euler: he is our master in everything. 

- Pierre Simon de Laplace 



Leonhard Euler was a Swiss mathematician who made signifi- 
cant contributions to a wide range of mathematics and physics 
including calculus and celestial mechanics (see [Eul-num] and 
[Eu2-num] for further details). 

The goal is to find an approximate solution to the problem 

y' = fix,y), yia)=c, (1.5) 

where f{x,y) is some given function. We shall try to approxi- 
mate the value of the solution x = b, where 6 > a is given. 
Sometimes such a method is called "numerically integrating 
(1.5)". 

Note: the first order DE must be in the form (1.5) or the 
method described below does not work. A version of Euler's 
method for systems of 1-st order DEs and higher order DEs will 
also be described below. 

Euler's method 

Geometric idea: The basic idea can be easily expressed in 
geometric terms. We know the solution, whatever it is, must go 
through the point (a, c) and we know, at that point, its slope is 



m = f{a,c). Using the point-slope form of a line, we conclude 
that the tangent line to the solution curve at (a, c) is (in {x, y)- 
coordinates, not to be confused with the dependent variable y 
and independent variable x of the DE) 

y = c-\- {x — a)f{a, c). 

In particular, if /z. > is a given small number (called the in- 
crement) then taking x = a-\-h the tangent-line approximation 
from calculus I gives us: 

y{a + h) = c-\-h ■ f{a, c). 

Now we know the solution passes through a point which is 
"nearly" equal to {a -\- h,c -\- h ■ f{a,c). We now repeat this 
tangent-line approximation with (a, c) replaced by (a + /i, c + 
h • f{a,c). Keep repeating this number-crunching at x = a, 
x = a-\-h, x = a-\- 2h, until you get to x = b. 

Algebraic idea: The basic idea can also be explained "alge- 
braically" . Recall from the definition of the derivative in calculus 
1 that 

_ y{x + h)- y{x) 

y{x) = , 

/i > is a given and small. This an the DE together give 
/(x, y{x)) = mNiMzM^) _ Now solve for y{x + h): 

y{x + h)= y{x) + h • f{x, y{x)). 

If we call h-f{x, y{x)) the "correction term" (for lack of anything 
better), call y{x) the "old value of ?/", and call y{x + h) the "new 
value of y" , then this approximation can be re-expressed 



Vnew = Void + h ■ f{x, yold)- 

Tabular idea: Let n > be an integer, which we call the 
step size. This is related to the increment by 

h — a 

h = . 

n 

This can be expressed simplest using a table. 



X 


y 


hf{x,y) 


a 


c 


hf{a,c) 


a-\- h 


c + hf{a, c) 




a + 2h 






b 


??? 


XXX 



The goal is to fill out all the blanks of the table but the xxx 
entry and find the ??? entry, which is the Euler's method 
approximation for y{b). 

Example 1.5.1. Use Euler's method with h = 1/2 to approxi- 
mate y{l), where 

y'-y = ^X-5, y{0) = l. 

Putting the DE into the form (1.5), we see that here f{x,y) = 
5x -\- y — 5, a = 0; c = 1. 



X 


y 


hf{x,y) = '-^ 







1 


-2 


1/2 


1 + 


(-2) = -1 


-7/4 


1 




-7/4) = -11/4 





so y{l) = — ^ = —2.75. This is the final answer. 

Aside: For your information, y = — 5x solves the DE and 
y(l) = e-5 = -2.28.... 

Here is one way to do this using SAGE ; 

SAGE 



sage : x, y=PolynomialRing (QQ, 2 , "xy " ) . gens ( ) 
sage: eulers_method {5*x+y-5, 1, 1, 1/3, 2) 

X y h*f (x, y) 

1 1 1/3 
4/3 4/3 1 
5/3 7/3 17/9 

2 38/9 83/27 
sage : eulers_method {5*x+y-5, 0, 1, 1/2, 1, method="none" ) 
[[0, 1], [1/2, -1], [1, -11/4], [3/2, -33/8]] 

sage: pts = eulers_method (5*x+y-5, 0, 1, 1/2, l,method="none"; 

sage: P = list_plot (pts) 

sage: show(P) 

sage: P = line (pts) 

sage: show(P) 

sage: PI = list_plot (pts ) 

sage: P2 = line (pts) 

sage: show(Pl+P2) 



The plot is given below. 




Figure 1.9: Euler's method with h = 1/2 for x' + x = 1, x{0) = 2. 



Improved Euler's method 

Geometric idea: The basic idea can be easily expressed in 
geometric terms. As in Euler's method, we know the solution 
must go through the point (a, c) and we know its slope there 
is m = /(a,c). If we went out one step using the tangent line 
approximation to the solution curve, the approximate slope to 
the tangent line at a: = a + h^y = c + h ■ f{a, c) would be 
m' = f{a+h, c+h-f{a, c)). The idea is that instead of using m = 
f{a, c) as the slope of the line to get our first approximation, use 
"^"t"^ . The "improved" tangent-line approximation at (a,c) is: 



V{a+h) - c+t^^ = ^^^ f{a,c)+f(a + h,c + h-f(a,c)) 



(This turns out to be a better apprpximation than the tangent- 
line approximation y{a + h) = c + h ■ f{a,c) used in Euler's 



method.) Now we know the solution passes through a point 
which is "nearly" equal to {a -\- h, c + h • ). We now repeat 
this tangent-line approximation with (a, c) replaced by (a+/i, c+ 
h • f{a,c). Keep repeating this number-crunching ai x = a, 
x = a-\-h, x = a-\- 2h, until you get to x = b. 

Tabular idea: The integer step size n > is related to the 
increment by 

b — a 

h = , 

n 

as before. 

The improved Euler method can be expressed simplest using a 
table. 



X 


y 


I^m+m' _ i^f{x,y)+f{x+h,y+h-f{x,y)) 


a 


c 


^ f{a,c}+f(a+h,c+h-f{a,c}) 


a-\- h 






a + 2h 






b 


??? 


XXX 



The goal is to fill out all the blanks of the table but the xxx 
entry and find the ??? entry, which is the improved Euler's 
method approximation for y{b). 

Example 1.5.2. Use the improved Euler's method with h = 1/2 
to approximate y{l), where 



y' - y = 5x - y{0) = 1. 



Putting the DE into the form (1-5), we see that here f{x,y) = 
5x -\- y — 5, a = 0; c = 1. We first compute the "correction 
term": 



h 



f{x,y)+f{x+h,y+h-f{x,y)) 
2 



5x+y-5+5{x+h)+{y+h-f{x,y))-5 
4 

5x+y-5+5{x+h)+{y+h-{5x+y-5)-5 



(l + |)5x + (l + |)?/ 
25x/4 + 5y/4-5. 



X 


y 


um+m' 25x+5y-10 
2 ~ 4 




1/2 
1 


-7/8 + (- 


1 

-15/8) = - 
-95/64) = 


-7/8 
-151/64 


-15/8 
-95/64 



so y{l) = — ^ = —2.35... This is the final answer. 

Aside: For your information, this is closer to the exact value 
y{\) = e — 5 = —2.28... than the "usual" Euler's method approx- 
imation of —2.75 we obtained above. 



Euler's method for systems and higher order DEs 



We only sketch the idea in some simple cases. Consider the 
DE 



y" + P{x)y' + q{x)y = fix), y{a) = ei, y'{a) = 62, 
and the system 

Vi = fi{x,yi,y2), yi{a) = ci, 

y2 = f2{x,yuy2), y2{a) = C2. 



We can treat both cases after first rewriting the DE as a system: 
create new variables yi = y and let y2 = y' ■ It is easy to see that 

?/'i = ?/2, ?/i(a) = ei, 

y'2 = fix) - q{x)yi - pix)y2, y2[a) = e^. 

Tabular idea: Let n > be an integer, which we call the 
step size. This is related to the increment by 

h — a 



This can be expressed simplest using a table. 



X 


yi 


hh{x,yi,y2) 


y2 


hf2{x,yuy2) 


a 


ei 


/l/l(a,6i,62) 




/i/2(a,6i,62) 


a-\- h 


ei + /i/i(a, 61,62) 




61 + /i/l(a,6i,62) 




a-\-2h 










b 


??? 


XXX 


XXX 


XXX 



The goal is to fill out all the blanks of the table but the xxx 
entry and find the ??? entries, which is the Euler's method 
approximation for y{b). 

Example 1.5.3. Using 3 steps of Euler's method, estimate x{\) , 
where x" - 2,x' + 2x = l, x(0) = 0, x'(0) = 1 

First, we rewrite x" — ?>x' + 2x = 1, x{Q) = x'{Q) = \, as a 
system of V* order DEs with ICs. Let xi = x, X2 = x' , so 



x'l = X2, Xi{0) = 0, 

X2 = 1 — 2xi + 3x2, ^2(0) = 1. 



This is the DE rewritten as a system in standard form. (In 
general^ the tabular method applies to any system but it must be 
in standard form.) 

Taking h = {1 — 0)/3 = 1/3, we have 



t 




X2/3 


X2 


(1 


— 2x1 + 3x2) /3 








1/3 


1 




4/3 


1/3 


1/3 


7/9 


7/3 




22/9 


2/3 


10/9 


43/27 


43/9 




XXX 


1 


73/27 


tXjtXjJu 






XXX 



So x{l) = xi{l) ~ 73/27 = 2.7.... 

Here is one way to do this using SAGE ; 
SAGE 

sage: RR = RealField (sci_not=0, prec=4, rnd='RNDU') 
sage: t, x, y = PolynomialRing (RR, 3, "txy" ) . gens ( ) 
sage: f = y; g = l-2*x+3*y 

sage: L = eulers_method_2x2 ( f , g, , , 1 , 1/ 3, 1 , method="none" ) 
sage: L 

[[0, 0, 1], [1/3, 0.35, 2.5], [2/3, 1.3, 5.5], 

[1, 3.3, 12], [4/3, 8.0, 24]] 
sage: eulers_method_2x2 (f , g, 0, 0, 1, 1/3, 1) 



t 


X 




h*f (t,x,y) 


y 


h*g(t,x,y) 










. 35 


1 


1 . 4 


1/3 





35 


.88 


2.5 


2 . 8 


2/3 


1 


3 


2.0 


5.5 


6.5 


1 


3 


3 


4.5 


12 


11 



sage: PI = list_plot { [ [p [ ] , p [ 1 ] ] for p in L] ) 
sage: P2 = line ( [ [p [ ] , p [ 1 ] ] for p in L] ) 
sage: show(Pl+P2) 



The plot of the approximation to x{t) is given below. 




Figure 1.10: Euler's method with h = 1/3 for x" — 3x' + 2x = 1, x{0) = 0, 
a;'(0) = 1. 

Exercise: Use SAG Band Euler's method with /i = 1/3 for the 
following problems: 
(a) Find the approximate values of x{l) and y{l) where 




x{0) = 1. 



1.6 Newtonian mechanics 



We briefly recall how the physics of the falling body problem 
leads naturally to a differential equation (this was already men- 
tioned in the introduction and forms a part of Newtonian me- 
chanics [M-mech]). Consider a mass m falling due to gravity. 
We orient coordinates to that downward is positive. Let x{t) 
denote the distance the mass has fallen at time t and v{t) its 
velocity at time t. We assume only two forces act: the force due 
to gravity, Fgrav^ and the force due to air resistence, F^es- In 
other words, we assume that the total force is given by 

-^total ~ Fgrav ~H Fres- 

We know that Fg^av = fng^ where p > is the gravitational 
constant, from high school physics. We assume, as is common 
in physics, that air resistance is proportional to velocity: F^es = 
—kv = —kx'{t), where A; > is a constant. Newton's second 
law [N-mech] tells us that Ftotai = "^cl = mx"{t). Putting these 
all together gives mx"{t) = mg — kx'{t), or 

v'it) + -v{t)=g. (1.6) 

m 

This is the differential equation governing the motion of a falling 
body. Equation (1.6) can be solved by various methods: separa- 
tion of variables or by integrating factors. If we assume v{0) = vq 
is given and if we assume A; > then the solution is 

„(t) = 7 + («o - ^)e-"/'». (1.7) 
In particular, we see that the limiting velocity is vumit = 



Example 1.6.1. Wile E. Coyote (see [W-mechJ if you haven't 
seen him before) has mass 100 kgs (with chute). The chute is 
released 30 seconds after the jump from a height of 2000 m. The 
force due to air resistence is given by F^es = —kv, where 



k = 

Find 



15, chute closed, 
100, chute open. 



(a) the distance and velocity functions during the time when 
the chute is closed (i.e., < t < 30 seconds), 

(b) the distance and velocity functions during the time when 
the chute is open (i.e., 30 < t seconds), 

(c) the time of landing, 

(d) the velocity of landing. (Does Wile E. Coyote survive the 
impact?) 

soln; Taking m = 100, g = 9.8, k = lb and v{0) = in (1.7), 
we find 

, , 196 196 

yi{t) = ^ 20*. 

This is the velocity with the time t starting the moment the 
parachutist jumps. After t = 30 seconds, this reaches the ve- 
locity vq = ^ — ^e~^/^ = 64.607.... The distance fallen is 

xi{t) = JqVi{u) du 

_ 196 ^ I 3920 ^-^t 3920 
- + 9~- 

After 30 seconds, it has fallen Xi{30) = ^ + Sg-^/^ = 
1529.283... meters. 



Taking m = 100, g = 9.8, k = 100 and v{0) = vq, we find 
, , 49 , /833 196 q/oA 

This is the velocity with the time t starting the moment Wile 
E. Coyote opens his chute (i.e., 30 seconds after jumping) . The 
distance fallen is 

X2[t) = Jq V2{u) du + Xi{30) 

_ 49 J. _ 833 -t I 196 p-ip-9/2 _|_ 71099 _|_ 3332 p-9/2 
5 ^ 15^ '3^^ ~'~45~'~9^ 



Now let us solve this using SAGE . 

SAGE - 



sage: RR = RealField (sci_not=0, prec=50, rnd='RNDU') 

sage : t = var ( ' t' ) 

sage: v = f unction {' v' , t) 

sage: m = 100; g = 98/10; k = 15 

sage: de = lambda v: m*diff (v,t) + k;*v - m*g 

sage: desolve_laplace (de (v(t) ) , ["t", "v"] , [0,0] ) 

'196/3-196*%e'-{3*t/20)/3' 

sage: solnl = lambda t: 196/3-196*exp (-3*t/20) /3 
sage: PI = plot (solnl (t) , 0, 30, plot_points=1000) 
sage: RR (solnl (30) ) 
64 . 607545559502 



This solves for the velocity before the coyote 's chute is opened, 
< t < 30. The last number is the velocity Wile E. Coyote is 
traveling at the moment he opens his chute. 

SAGE 

sage : t = var (' t' ) 

sage: v = function ('v', t) 

sage: m = 100; g = 98/10; k = 100 

sage: de = lambda v: m*diff (v,t) + k*v - m*g 

sage: desolve_laplace (de (v(t) ) , ["t", "v"] , [0,RR(solnl (30) ) ] ) 

' 631 931 *%e'^-t/ 11530 + 4 9/5' 

sage: soln2 = lambda t: 49/5+ (631931/11530) *exp (- (t-30) ) 

+ solnl(30) - (631931/11530) - 49/5 
sage: RR { soln2 (30 ) ) 
64 . 607545559502 
sage: RR ( solnl (30 ) ) 
64 . 607545559502 

sage: P2 = plot {soln2 (t) , 30, 50,plot_points=1000) 
sage: show(Pl+P2) 



This solves for the velocity after the coyote 's chute is opened, t > 



30. The last command plots the velocity functions together as a 
single plot. (You would see a break in the graph if you omitted 
the SAGE '5 plot option ,plot_points=1000. That is because 
the number of samples taken of the function by default is not 
sufficient to capture the jump the function takes att = 30.) The 
terms at the end o/soln2 were added to insure X2{S0) = xi{30). 
Next, we find the distance traveled at time t: 

SAGE 

age: integral { solnl (t ), t ) 
3920*e' (- (3*t/20) ) /9 + 196*t/3 

sage: xl = lambda t: 3920*e' (- (3*t/20) ) /9 + 196*t/3 
sage: RR(xl(30)) 
1964 . 8385851589 



This solves for the distance the coyote traveled before the chute 
was open, < t < 30. The last number says that he has gone 
about 1965 meters when he opens his chute. 
SAGE 

sage: integral (soln2 (t) , t) 

49*t/5 - {631931*e' (30 - t)/11530) 

sage: x2 = lambda t: 49*t/5 - (631931*e" (30 - t)/11530) 

+ xl(30) + (631931/11530) - 49*30/5 
sage: RR(x2(30.7)) 
1999.2895090436 
sage: P4 = plot (x2 (t) , 30, 50) 
sage: show(P3+P4) 



(Again, you see a break in the graph because of the round-off 
error.) The terms at the end ofx2 were added to insure X2(30) = 
Xi{30). You know he is close to the ground att = 30; and going 



quite fast (about 65 m/sl). It makes sense that he will hit the 
ground soon afterwards (with a large puff of smoke, if you've 
seen the cartoons), even though his chute will have slowed him 
down somewhat. 

The graph of the velocity < t < 50 is in Figure 1.11. Notice 
how it drops at t = 30 when the chute is opened. The graph of 
the distance fallen < t < 50 is in Figure 1.12. Notice how it 
slows down at t = 30 when the chute is opened. 




Figure 1.11: Velocity of falling parachutist. 



The time of impact is timpact = 30.7.... This was found numer- 
ically by a "trial- and- error" method of solving X2{t) = 2000. 
The velocity of impact is V2{timpact) ~ 37 m/s. 

Exercise: Drop an object with mass 10 kgs from a height of 
2000 m. Suppose the force due to air resistence is given by 
Fres = —^Ov. Find the velocity after 10 seconds using SAGE. 
Plot this velocity function for < t < 10. 



Figure 1.12: Distance fallen by a parachutist. 

1.7 Application to mixing problems 

Suppose that we have two chemical substances where one is 
soluable in the other, such as salt and water. Suppose that we 
have a tank containing a mixture of these substances, and the 
mixture of them is poured in and the resulting "well-mixed" 
solution pours out through a value at the bottom. (The term 
"well-mixed" is used to indicate that the fluid being poured in 
is assumed to instantly dissolve into a homogeneous mixture the 
moment it goes into the tank.) The crude picture looks like this: 

Assume for concreteness that the chemical substances are salt 
and water. Let 

• A{t) denote the amount of salt at time t, 

• FlowRateIn = the rate at which the solution pours into the 
tank. 




Figure 1.13: Solution pours into a tank, mixes with another type of solution, 
and then pours out. 



• FlowRateOut = the rate at which the mixture pours out of 
the tank, 

• Cin = "concentration in" = the concentration of salt in the 
solution being poured into the tank, 

• Cout = "concentration out" = the concentration of salt in 
the solution being poured out of the tank, 

• Rin = rate at which the salt is being poured into the tank 
= (FlowRateIn)(Cin), 

• Rout = rate at which the salt is being poured out of the 
tank = (FlowRateOut) (Cout)- 



Remark 1.7.1. Some things to make note of: 



• // FlowRateIn = FlowRateOut then the "water level" of the 
tank stays the same. 

• We can determine Cout cls a function of other quantities: 

A(t} 



Cout — 



Tit) 



where T{t) denotes the volume of solution in the tank at 
time t. 

The rate of change of the amount of salt in the tank, A'{t), 
more properly could he called the "net rate of change" . If 
you think if it this way then you see A'[t) = Rin — Rout- 

Now the differential equation for the amount of salt arises from 
the above equations: 

A'{t) = (FlowRateIn) - (FlowRateOut) 

1 [tj 

Example 1.7.1. Consider a tank with 200 liters of salt-water 
solution, 30 grams of which is salt. Pouring into the tank is a 
brine solution at a rate of 4 liters/minute and with a concentra- 
tion of 1 grams per liter. The "well-mixed" solution pours out 
at a rate of 5 liters /minute. Find the amount at time t. 
We know 



A'{t) = (FlowRateIn)Q„-(FlowRateOut)^ = 4-5 ^^^^^^ , A{0) 

Writing this in the standard form A' + pA = q, we have 

_ J iJ.(t)q{t) dt + C 



where ji = e^^^^^'^* = e~^^ 2ohi^t = (200 is the "integrating 
factor". This gives A{t) = 200 - t + C • (200 - ^)^ where the 
initial condition implies C = — 170 • 200~^. 

Here is one way to do this using SAGE ; 

SAGE 

sage : t = var { ' t' ) 

sage: A = f unction (' A' , t) 

sage: de = lambda A: diff{A, t) + (5/ (200-t) ) *A - 4 

sage : de solve (de ( A (t ) ) , [ A, t ] ) 
' {%c-l/ {t-200) ~4) * {t-200) ^5' 



This is the form of the general solution. (SAGE uses Maxima 
and %c is Maxima's notation for an arbitrary constant.) Let us 
now solve this general solution for c, using the initial conditions. 

SAGE 

sage : c = var { ' c' ) 

sage: solnA = lambda t: (c - 1/ (t-200) "4) * (t-200) "5 
sage: solnA(t) 

(c - (l/(t - 200)'4))*(t - 200)^5 

sage: solnA(O) 

-320000000000* (c - 1/1600000000) 
sage: solve ([ solnA ( ) == 30], c) 
[c == 17/32000000000] 
sage: c = 17/32000000000 
sage: solnA(t) 

(17/32000000000 - (1/ (t - 200)"4))*(t - 200)"5 
sage: P = plot (solnA (t) , 0, 200) 
sage: show(P) 



Figure 1.14: A{t), 0<t< 200, A' = 4 - 5A(t)/(200 - t), A{0) = 30. 



Exercise: Now use SAGE to solve the same problem but with 
the same flow rate out as 4 liters/min (so the "water level" in the 
tank is constant). Find and plot the solution A{t), < t < 200. 
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2.1 Linear differential equations 



We want to describe the form a solution to a linear ODE can 
take. Before doing this, we introduce two pieces of terminology. 

• Suppose fi{t), f2{t), . . . , fn{t) are given functions. A lin- 
ear combination of these functions is another fucntion of 
the form 

Cl/l(t)+C2/2(t) + ...,+C„/n(^), 

for some constants ci, c„. For example, 3 cos(t) — 2 sin(t) 
is a linear combination of cos(t), sin(t). 

• A linear ODE of the form 

+ h{t)y^^-'^ + ... + bn-iit)y' + bnit)y = /(t), (2.1) 

is called homogeneous if f{t) = (i.e., / is the function) 
and otherwise it is called non-homogeneous. 

The following result describes the general solution to a linear 
ODE. 

Theorem 2.1.1. Consider a linear ODE having of the form 
(2.1), for some given continuous functions bi{t), . . . , bn{t), and 
f{t). Then the following hold. 

• There are n functions yi{t), . . . , yn{t) (called fundamen- 
tal solutions^; each satisfying the homogeneous ODE 

y'"Uhi{t)yl-"-'> + ... + h„.iit)y' + h„(t)y = 0, l<i<n, 

(2.2) 



such that every solution to (2.2) is a linear combination of 
these functions yi, . . . ,yn. 



• Suppose you know a solution yp{t) (a particular solution^ 
to (2.1). Then every solution y = y{i) (the general solu- 
tion j to the DE (2.1) has the form 

y{t)=yh{t)^yp{t), (2.3) 

where yu (the "homogeneous part" of the general solution) 
is a linear combination 

Vhit) = ciyi{t) + y2it) + ... + Cnynit), 
for some constants Ci, I < i < n. 

• Conversely, every function of the form (2.3), for any con- 
stants Ci for 1 <i <n, is a solution to (2.1). 

Example 2.1.1. Recall the example in the introduction where 
we looked for functions solving x' -\- x = 1 by "guessing". The 
function Xp{t) = 1 is a particular solution to x' -\- x = 1. The 
function xi{t) = e~* is a fundamental solution to x' -\- x = 0. 
The general solution is therefore x{t) = l-\-cie~*, for a constant 
ci- 

Example 2.1.2. The charge on the capacitor of an RLC elec- 
trical circuit is modeled by a 2-nd order linear DE [C-linear]. 
Series RLC Circuit notations: 

• E = E{i) - the voltage of the power source (a battery or 
other "electromotive force" , measured in volts, V) 



• q = q{t) - the current in the circuit (measured in coulombs, 
G) 

• i = i{t) - the current in the circuit (measured in amperes, 
A) 

• L - the inductance of the inductor ( measured in henrys, H) 

• R - the resistance of the resistor ( measured in ohms, Q ); 

• C - the capacitance of the capacitor (measured in farads, 
F) 

The charge q on the capacitor satisfies the linear IPV: 
Lq"+Rq' + ^q=E(t), q(Q) = q», (/'(O) = io- 



Figure 2.1: RLC circuit. 

Example 2.1.3. Recall the example in the introduction where 
we looked for functions solving x' + x = 1 by "guessing" . The 
function Xp{t) = I is a particular solution to x' + x = 1. The 
function xi{t) = e^* is a fundamental solution to x' + x = 0. 
The general solution is therefore x{t) = l + cie~^ , for a constant 



Example 2.1.4. The displacement from equilibrium of a mass 
attached to a spring is modeled by a 2-nd order linear DE [0-ivpJ. 
SSpring-mass notations: 

• f{t) - the external force acting on the spring (if any) 

• X = x{t) - the displacement from equilibrium of a mass 
attached to a spring 

• m - the mass 

• b - the damping constant (if, say, the spring is immersed in 



a fluid) 



• k - the spring constant. 



The displacement x satisfies the linear IPV: 



mx" + bx' -\- kx = fit) 



x{fS) = xo, x'ifS) = Vq. 




Figure 2.2: spring-mass model. 



Notice that each general solution to an n-th order ODE has 
n "degrees of freedom" (the arbitrary constants q). According 
to this theorem, to find the general solution of a linear ODE, 
we need only find a particular solution yp and n fundamental 
solutions yi{t), . . . , 

Example 2.1.5. Let us try to solve 

x' -\- X = 1, x{0) = c, 

where c = 1, c = 2, and c = 3. (Three different IVP's, three 
different solutions, find each one.) 

The first problem, x' -\- x = 1 and x(0) = 1, is easy. The 
solutions to the DE x' -\- x = 1 which we "guessed at" in the 
previous example, x{t) = 1, satisfies this. 

The second problem, x' -\-x = 1 and x(0) = 2, is not so simple. 
To solve this (and the third problem), we really need to know 
what the form is of the "general solution". 

According to the theorem above, the general solution x has the 
form X = Xp-\- Xh- In this case, Xp(t) = 1 and Xh{t) = ciXi(t) = 
cie~^ , by an earlier example. Therefore, every solution to the 
DE above is of the form x{t) = 1 + cie~^ , for some constant ci. 
We use the initial condition to solve for c\ : 

• a;(0) = 1; 1 = x(0) = 1 + c\^ = 1 + ci so ci = and 
x{t) = 1. 

• x{0) = 2: 2 = x{0) = 1 + cie^ = 1 + ci so ci = 1 and 

x{t) = 1 + e-K 

• a:(0) =3:3 = x(0) = 1 + cie^ = 1 + ci so ci = 2 and 

x{t) = 1 + 2e-K 



Here is one way to use SAGE to solve for ci. (Of course, you 
can do this yourself but this shows you the SAGE syntax for 
solving equations. Type solve? in SAGE to get more details.) 
We use SAGE to solve the last IVP discussed above and then to 
plot the solution. 

SAGE 

sage : t = var ( ' t' ) 
sage: cl = var('cl') 

sage: solnx = lambda t: l+cl*exp(-t) 
sage: solnx(O) 
cl + 1 

sage: solve ([ solnx ( ) == 3],cl) 

[cl == 2] 

sage : c_l = solve ( [solnx (0) == 3],cl)[0]. rhs ( ) 

sage: c_l 

2 

sage: solnxl = lambda t: l+c*exp(-t) 
sage: plot (solnxl (t) , 0, 2) 

Graphics object consisting of 1 graphics primitive 
sage: P = plot (solnxl (t) , 0, 2 ) 
sage: show{P) 

sage: P = plot (solnxl (t) , 0, 5) 
sage: show(P) 



This plot is shown below. 



1 2 3 A B 



Figure 2.3: Solution to IVP x' + x = 1, x(0) = 3. 



Exercise; Use SAGE to solve and plot the solution to x' + x = 1 
and x{0) = 2. 



2.2 Linear differential equations, continued 



To better describe the form a solution to a linear ODE can 
take, we need to better understand the nature of fundamental 
solutions and particular solutions. 
Recall that the general solution to 

2/^"^ + h{t)y^^-'^ + ... + bn-i{t)y' + br,{t)y = f{t), 

has the form y = yp + yn^ where yn is a linear combination of 
fundamental solutions. For example, the general solution to the 
spring-mass equation 

x" -\-x = l 

has the form x = x{t) = 1 -\- ci cos(t) + C2 sin(t) for the displace- 
ment from the equilibrium position. Suppose we are also given n 
initial conditions y{xQ) = ao, y'ixo) = ai, . . . , y^'^~^\xQ) = a„_i. 
For example, we could impose the initial position and initial 
velocity on the mass: x{0) = xq and x'{0) = vq. Of course, 
no matter what xq and vq are are given, we want to be able to 
solve for the coefficients ci, C2 in x{t) = 1 + ci cos(t) + C2 sin(t) 
to obtain a unique solution. More generally, we want to be able 
to solve an n-th order IVP and obtain a unique solution. A few 
questions arise. 

• How do we know this can be done? 

• How do we know that there isn't a linear combination of 
fundamental solutions which isn't (i.e., the zero function)? 

The complete answer actually involves methods from linear 
algebra which go beyond this course. The basic idea though 



is not hard to understand and it involves what is cahed "the 
Wronskian^" [W-hnear]. We'll have to explain what this means 
first. If f2{t), fn{t) are given n-times differentiable 

functions then their fundamental matrix is the matrix 



/ flit) f2{t) 



fn{t) 

m 



\ 



The determinant of the fundamental matrix is called the Wron- 
skian, denoted W{J\, ...,/«)• The Wronskian actually helps us 
answer both questions above simultaneously. 

Example 2.2.1. Take fi{t) = sin^{t), f2{t) = cos^(t), and 
f^{t) = 1. SAGE allows us to easily compute the Wronskian: 

SAGE 



sage: SR = SymbolicExpressionRing ( ) 

sage: MS = MatrixSpace (SR, 3, 3 ) 

sage: Phi = MS { [ [ sin (t ) " 2 , cos (t ) " 2 , 1 ] , 

[dif f (sin (t) "2, t) , dif f (cos (t) "2, t) , 0] , 
[diff (sin (t) '•2,t,t) ,diff (cos (t)''2,t,t),0]] 

sage: Phi 



[ sin(t)"2 cos(t)"2 

[ 2*cos (t) *sin (t) -2*cos (t) *sin (t) 

[2*cos(t)'2 - 2*sin(t)'2 2*sin(t)'2 - 2*cos (t) '2 
sage : Phi . det ( ) 




1] 
0] 
0] 



^ Josef Wronski was a Polish-born French mathemtician who worked in many different 
areas of appUed mathematics and mechanical engineering [Wr-linear]. 



Here Phi.detO is the determinant of the fundamental matrix 
Phi. Since it is zero, this means W {sin(ty , cos{tY , 1) = 0. 
(Note: the above entry for Phi should all be on one line. For 
typographical reasons, we have spread it out to 3 lines.) The 
entries for the symbolic expression ring SR and the 3x3 ma- 
trix space MS above are used to construct the matrix Phi having 
symbolic entries. 



We try one more example: 

SAGE 



sage: SR = SymbolicExpressionRing ( ) 

sage: MS = MatrixSpace (SR, 2, 2 ) 

sage: Phi = MS ( [ [ sin (t ) " 2 , cos (t ) " 2 ] , 

[diff (sin (t) "2, t) , dif f (cos (t) "2, t) ] ] ) 

sage: Phi 



[ sin (t) "^2 cos (t) "2] 

[ 2*cos (t) *sin (t) -2*cos (t) *sin (t) ] 
sage : Phi . det ( ) 

-2*cos (t) *sin (t) '3 - 2*cos (t) '3*sin (t) 



This means W{sin{ty, COS (t)'^) = — 2cos(t) sin(t)^— 2 cos(t)^ sin(i), 
which is non-zero. 

If there are constants Ci, c„, not all zero, for which 



cifiit) + C2f2it)--- + Cnfnit) = 0, for alH, (2.4) 

then the functions fi {I < i < n) are called linearly depen- 
dent. If the functions fi (I < i <n) are not linearly dependent 
then they are called linearly independent (this definition is 
frequently seen for linearly independent vectors [L-linear] but 
holds for functions as well). This condition (2.4) can be inter- 
preted geometrically as follows. Just as cix + C2y = is a line 
through the origin in the plane and Cix + C2y + Csz = is a plane 
containing the origin in 3-space, the equation 

CiXi + C2X2 h CnXn = 0, 

is a "hyperplane" containing the origin in n-space with coordi- 
nates {xi,...,Xn)- This condition (2.4) says geometrically that 



the graph of the space curve f{t) = . . . , fn{t)) lies en- 

tirely in this hyperplane. If you pick n functions "at random" 
then they are "probably" linearly independent, because "ran- 
dom" space curves don't lie in a hyperplane. But certainly not 
all collections of functions are linearly independent. 

Example 2.2.2. Consider just the two functions fi{t) = sin^(t), 
f2{t) = cos^(t). We know from the SAGE computation in the 
example above that these functions are linearly independent. 

SAGE 

sage: P = parametric_plot { (sin (t) "2, cos (t) "2) , 0, 5) 
sage: show(P) 



The SAGE plot of this space curve r{t) = (sin(t)^, cos(t)^) is 
given below. It is obviously not contained in a line through the 
origin, therefore making it geometrically clear that these func- 
tions are linearly independent. 

The following two results answer the above questions. 

Theorem 2.2.1. fWronskian test; /2(t); fn{t) 

are given n-times differentiable functions with a non-zero Wron- 
skian then they are linearly independent. 

As a consequence of this theorem, and the SAGE computation 
in the example above, fi(t) = sin^(i), /2(t) = cos^(t), are lin- 
early independent. 

Theorem 2.2.2. Given any homogeneous n-th linear ODE 

+ hi{t)y^^-'^ + ... + hn-i{t)y' + hn{t)y = 0, 

with differentiable coefficients, there always exists n solutions 
y\{t), yn{t) which have a non-zero Wronskian. 



Figure 2.4: Parametric plot of (sin(i(:)^, cos(t)^). 

The functions yi{t), yn{t) in the above theorem are called 
fundamental solutions. 

We shall not prove either of these theorems here. Please see 
[BD-intro] for further details. 

Exercise: Use SAGE to compute the Wronskian of 

(a) /i(t) = sin(t), /2(t)=cos(^), 

(b) m = 1, hit) = t, hit) = t\ hit) = t'. 
Check that 

(a) yiit) = sin(t), 2/2 (^) = cos(i) are fundamental solutions for 

y" + y = 0, 

(d) yiit) = 1, y2it) = t, ysit) = t^, yiit) = t^ are fundamental 
solutions for y'^'^^ = y"" = 0. 



2.3 Undetermined coefficients method 



The method of undetermined coefficients [U-uc] can be used to 
solve the following type of problem. 

PROBLEM: Solve 

ay" + hy' ^cy = f{x), (2.5) 

where a 7^ 0, 6 and c are constants and x is the independent 
variable. (Even the case a = can be handled similarly, though 
some of the discussion below might need to be slightly modified.) 
Where we must assume that f{x) is of a special form. 

More-or-less equivalent is the method of annihilating operators 
[A-uc] (they solve the same class of DEs), but that method will 
be discussed separately. 

For the moment, let us assume f{x) has the form ai ■ p{x) • 
Qa2x . cos(a3x), or ai ■ p{x) ■ e"^^ • sin(a3x), where ai, a2, are 
constants and p{x) is a polynomial. 

Solution: 

• Find the "homogeneous solution" y^ to ay" + by' + cy = 
0, yh = ciyi + C2?/2- Here yi and y2 are determined as 
follows: let ri and r2 denote the roots of the characteristic 
polynomial aD^ + bD + c = 0. 

— ri 7^ r2 real: set yi = e'"^^, y2 = e^'^^. 

— ri = r2 real: if r = ri = r2 then set yi = e*"^, y2 = xe^^. 

— ri, r2 complex: if ri = a + i/?, r2 = a — i(3^ where a and 
(3 are real, then set yi = e"^ cos (/3a;), y2 = e"^ sin(/3x). 



• Compute f{x), f'{x), f"{x), ... . Write down the list of 
all the different terms which arise (via the product rule), 
ignoring constant factors, plus signs, and minus signs: 

If any one of these agrees with yi or y2 then multiply them 
all by X. (If, after this, any of them still agrees with yi or 
y2 then multiply them all again hy x.) 

• Let yp be a linear combination of these functions (your 
"guess"): 

Vp = Aiti{x) + ... + Aktkix). 

This is called the general form of the particular solu- 
tion. The Aj's are called undetermined coefficients. 

• Plug yp into (2.5) and solve for Ai, A^. 

• hei y = yu + yp = yp + cm + 022/2- This is the general 
solution to (2.5). If there are any initial conditions for 
(2.5), solve for then ci,C2 now. 



Diagramatically: 



Factor characteristic polynomial 
i 



Compute yfi 



Compute the general form of the particular, 

i 

Compute the undetermined coefficients 



Answer: y = yf^-\- y^. 



Examples 

Example 2.3.1. Solve 

y" — y = cos(2x). 

• The characteristic polynomial is r^ — 1 = Q, which has ±1 for roots. The 
"homogeneous solution" is therefore y^ = cie^ + C2e~^. 

• We compute f{x) = cos{2x), f'{x) = — 2sin(2x), f"{x) = — 4cos(2a;), ... . 
They are all linear combinations of 

fi{x) = cos(2a;), /2(x) = sin(2a;). 
None of these agrees with yi = e^ or y2 = e~^, so we do not multiply by x. 

• Let yp be a linear combination of these functions: 

yp = Ai cos(2x) + A2 sin(2x). 

• You can compute both sides of y'p — yp = cos(2x); 

(-4^1 cos(2a;) - 4.42 sin(2x)) - {Ai cos(2x) + A2 sin(2x)) = cos(2x). 

Equating the coefficients of cos{2x), sin(2x) on both sides gives 2 equations 
in 2 unknowns: —bAi = 1 and — 5^2 = 0. Solving, we get Ai = —1/5 and 

^2 = 0. 

• The general solution: y = yn + yp = cie^ + C2e^^ — | cos(2a;). 
Example 2.3.2. Solve 

y" — y = xcos(2x). 

• The characteristic polynomial is — 1 = 0, which has ±1 for roots. The 
"homogeneous solution" is therefore yh = cie^ + 026"^. 

• We compute f{x) = xcos(2x), f'{x) = cos(2x)— 2x sin(2x), f"{x) = — 2sm(2x)— 
2sin(2x) — 2xcos(2a;), ... . They are all linear combinations of 

fi{x) = cos(2a;), /2(a;) = sin(2x), fsix) = a:cos(2x), ./4(x) = a;sin(2a:). 
None of these agrees with yi = e^ or y2 = e~^ , so we do not multiply by x. 



• Let Up he a linear combination of these functions: 



yp = Ai cos(2x) + A2 sm(2a;) + ^3a;cos(2a;) + ^4a;sin(2x). 

• In principle, you can compute both sides ofy'^~yp = a:;cos(2a:) and solve for 
the Ai's. (Equate coeffi,cients o/.tcos(2.x) on both sides, equate coefficients 
o/cos(2x) on both sides, equate coefficients of xsin{2x) on both sides, and 
equate coefficients of sin(2a;) on both sides. This gives 4 equations in 4 
unknowns, which can be solved.) You will not be asked to solve for the Ai 's 
for a problem this hard. 

Example 2.3.3. Solve 

y" + 4y = a;cos(2a;). 

• The characteristic polynomial is r"^ + A = 0, which has ±2z for roots. The 
"homogeneous solution" is therefore y^ = c\ cos(2x) + C2 sin(2a;). 

• We compute fix) = xcos(2a;), f'[x) = cos(2x)— 2x siii(2a;), fix) = — 2sm(2x) — 
2sin(2a;) — 2xcos(2x), ... . They are all linear combinations of 

/i(x) = cos(2x), f2ix) = sin(2x), faix) = xcos(2x), -fiix) = xsm(2x). 

Two of these agree with yi = cos(2x) or 1/2 = sin(2x), so we do multiply by 
x: 

fiix) = a;cos(2x), f2ix) = a;sin(2a;), fsix) = a;^cos(2x), ./4(a;) = x^sm(2a;). 

• Let yp be a linear combination of these functions: 

yp = ^ia;cos(2a:) + ^22: sin(2a;) + Asx"^ cos(2a;) + A^x^ sin(2x). 

• In principle, you can compute both sides of y'^ + 4yp = x cos(2a;) and solve 
for the Ai 's. You will not be asked to solve for the Ai 's for a problem this 
hard. 



More generally, suppose that you want to solve ay" + by' + cy = 
/(x), where f{x) is a sum of functions of the above form. In 
other words, f{x) = fi{x) + f2ix) + ... + fkix), where each fj{x) 



is of the form c-p{x) • • cos{bx), or c-p{x) ■ e"^ • sm{bx), where 
a, b, c are constants and p{x) is a polynomial. You can proceed 
in either one of the following ways. 

1. Split up the problem by solving each of the k problems 
ay" + by' + cy = fj{x), 1 < j < k, obtaining the solution 
y = yj{x), say. The solution to ay" + by' + cy = f{x) is 
then y = yi -\- y2 -\- .. -\- yk (the superposition principle) . 

2. Proceed as in the examples above but with the following 
slight revision: 

• Find the "homogeneous solution" y^ to ay" + by' = 

cy = 0,yh = ciyi + C2I/2- 

• Compute f{x), f'{x), f"{x), ... . Write down the list 
of all the different terms which arise, ignoring constant 
factors, plus signs, and minus signs: 

t2ix), tkix). 

• Group these terms into their families. Each family is 
determined from its parent (s) - which are the terms in 
f{x) = fi{x) + f2{x) + ... + fk{x) which they arose form 
by differentiation. For example, if f{x) = xcos(2x) + 

sin(x) + sin(2x) then the terms you get from differ- 
entiating and ignoring constants, plus signs and minus 
signs, are 

X cos(2a;) , x sin(2a;) , cos(2a;) , sin(2a;) , (from x cos(2a;)) , 
sin(a;) , cos (x) , (from sin(a;) ) , 



and 

sin(2a;), cos(2a;), (from sin(2a:)). 

The first group absorbes the last group, since you can 

only count the different terms. Therefore, there are 

only two families in this example: {x cos(2x) , x sin(2a:) , cos(2a:) , sin(2rc) } 

is a "family" (with "parent" xcos(2x) and the other 

terms as its "children") and {e~^ sin(x), cos(a;)} is 

a "family" (with "parent" sin(x) and the other term 

as its "child"). 

If any one of these terms agrees with yi or y2 then 
multiply the entire family hj x. In other words, if any 
child or parent is "bad" then the entire family is "bad" . 
(If, after this, any of them still agrees with yi or y2 then 
multiply them all again by x.) 

• Let yp be a linear combination of these functions (your 
"guess"): 

Up = Aiti{x) + ... + Aktk{x). 

This is called the general form of the particular 
solution. The Aj's are called undetermined coeffi- 
cients. 

• Plug yp into (2.5) and solve for Ai, A^. 

• Let y = yh + yp = yp + ciyi + 022/2- This is the general 
solution to (2.5). If there are any initial conditions for 
(2.5), solve for then ci,C2 last - after the undetermined 
coefficients. 

Example 2.3.4. Solve 

y'" -y"-y' + y = i2xe^ 



We use SAGE /or this. 

SAGE 



sage: x = var("x") 

sage: y = f unction ( "y" , x) 

sage: R.<D> = PolynomialRing (QQ, "D") 

sage: f = D"3 - D"2 - D + 1 

sage : f . factor ( ) 

(D + 1) * (D - 1) '2 
sage : f . roots { ) 

[ (-1, 1) , (1, 2) ] 



5*0 the roots of the characteristic polynomial are 1.1,-1, which 
means that the homogeneous part of the solution is 

Vh = cie^ + C2a;e^ + cse""". 
SAGE 



sage: de = lambda y: diff(y,x,3) - diff{y,x,2) - diff(y,x,l) + y 
sage: cl = var("cl"); c2 = var("c2"); c3 = var("c3") 
sage: yh = cl*e'x + c2*x*e'x + c3*e' (-x) 
sage : de (yh) 


sage: de (x'3*e'x- (3/2) *x'2*e'x) 
12*x*e'x 



This just confirmed that yh solves y'" — y" — y' -\- 1 = 0. Using 
the derivatives of F{x) = 12xe^, we generate the general form 
of the particular: 

SAGE 

sage: F = 12*x*e"x 



sage: diff{F,x,l); diff(F,x,2); diff(F,x,3) 

12*x*e'x + 12*e'x 

12*x*e"x + 24*e"x 

12*x*e'x + 35*e'x 
sage: Al = var{"Al"); A2 = var("A2") 
sage: yp = Al*x"2*e'x + A2*x'3*e'x 



Now plug this into the DE and compare coefficients of like terms 
to solve for the undertermined coefficients: 

SAGE 



sage : de (yp) 

12*x*e'~x*A2 + 6*e'x*A2 + 4*e'x*Al 
sage: solve {[12*A2 == 12, 6*A2+4*A1 0],A1,A2) 

[ [Al == -3/2, A2 == 1] ] 



Finally, lets check if this is correct: 

SAGE 

sage: y = yh + (-3/2) *x"2*e"x + (1) *x"3*e'x 
sage : de (y) 
12*x*e'~x 



Exercise: Using SAGE , solve 

y'" -y" + y'-y = 12xe\ 

(Hint: You may need to replace sage : R.<D> = PolynomialRing(QQ, 
by sage: R.<D> = PolynomialRingCCC, "D").) 



2.3.1 Annihilator method 



PROBLEM: Solve 

ay" + hy' + cy = f{x). (2.6) 

We assume that f{x) is of the form c • p{x) ■ e^^ ■ cos(6a:), or 
c • p{x) ■ e"^ • sm{bx), where a, 6, c are constants and p{x) is a 
polynomial. 

soln: 

• Write the ODE in symbolic form {aD^ + bD + c)y = f{x). 

• Find the "homogeneous solution" y^ to ay" -\-by' = cy = 0, 
Vh = ciyi + C2?/2. 

• Find the differential operator L which annihilates f{x): 
Lf{x) = 0. The following table may help. 



fiiiictioii 


aiiiiiliilator 




jjk+i 




{D - a)^+i 


Qr.k^ax cos(/5x) 









• Find the general solution to the homogeneous ODE, L ■ 
[aD'^ + bD + c)y = 0. 

• Let yp be the function you get by taking the solution you 
just found and subtracting from it any terms in y^. 

• Solve for the undetermined coefficients in yp as in the method 
of undetermined coefficients. 



Example 

Example 2.3.5. Solve 

y" — y = cos(2x). 

• The DE is (D^ _ i^^^ = cos(2x). 

• The characteristic polynomial is — 1 = 0, which has ±1 for roots. The 
"homogeneous solution" is therefore y^ = cie^ + C2e~^. 

• We find L = D'^ + 4 annihilates cos(2x). 

• We solve {D"^ + 4)(D^ — l)y = 0. The roots of the characteristic polynomial 
(r^ + 4)(r^ — 1) are ±2i, ±1. The solution is 

y = Ai cos(2x) + A2 sm(2x) + Ase"" + A^e'"^ . 

• This solution agrees with in the last two terms, so we guess 

yp = Ai cos(2x) + A2 sin(2a;). 

• Now solve for Ai and A2 as before: Compute both sides ofy'p — yp = cos(2x), 

(-4^1 cos(2x) - 4^2 sm(2x)) - (^1 cos(2x) + A2 siii(2x)) = cos(2x). 

Next, equate the coefficients o/cos(2.t), sin(2a;) on both sides to get 2 equa- 
tions in 2 unknowns. Solving, we get Ai = —1/5 a,nd = 0. 

• The general solution: y = yn + yp = cie^ + C2e~^ — g cos(2x). 



2.4 Variation of parameters 



Consider an ordinary constant coefficient non- homogeneous 2nd 
order linear differential equation, 

ay" + hy' + cy = F{x) 

where F{x) is a given function and a, 6, and c are constants. (For 
the method below, a, 6, and c may be allowed to depend on the 
independent variable x as well.) Let yi{x)^ y2{x) be fundamental 
solutions of the corresponding homogeneous equation 

ay" + hy' + cy = 0. 

The method of variation of parameters is originally attributed 
to Joseph Louis Lagrange (1736-1813) [L-var]. It starts by as- 
suming that there is a particular solution in the form 

yp{x) = ui{x)yi{x) + U2{x)y2{x), 

where Ui{x), U2{x) are unknown functions [V-var]. 
In general, the product rule gives 

if 9)' = f'g + fg', 

ifg)" = f"g + 2 f'g' + 

if 9)"' = f'g + ^fg' + 3/y + //', 

and so on, following Pascal's triangle. 



1 

1 1 



and so on. 



Using SAGE , this can be check as follows: 

SAGE 



sage : t = var ( ' t' ) 

sage: x = function (' x' , t) 
sage: y = f unction (' y' , t) 
sage: dif f (x (t) *y (t) , t) 

x(t) *diff (y{t) , t, 1) + y (t) *diff (x(t) , t, 1) 
sage: dif f (x (t) *y (t) , t, t) 

X (t) *diff (y (t) , t, 2) + 2*diff(x(t), t, 1) *dif f (y (t) , t, 1) 

+ y (t) *dif f (x (t) , t, 2) 
sage: dif f (x (t) *y (t) , t, t, t) 

x(t) *diff (y (t) , t, 3) + 3*diff(x(t), t, 1) *dif f (y (t) , t, 2) 
+ 3*diff(x(t), t, 2) *diff (y (t) , t, 1) + y (t ) *dif f (x (t ) , t, 3) 



By assumption, yp solves the ODE, so 



After some algebra, this becomes: 

a{u\yi + U2I/2)' + a{u\y[ + u^y'i) + h{u[yi + ^21/2) = F. 
If we assume 



Cramer's rule says that the solution to this system is 



+ ^y'p + = F{x). 



u'lVi + W2?/2 = 
then we get massive simplification: 



a{u[y[ + M2?/2) = F. 




Note that the Wronskian of the fundamental solutions ^2) 
is in the denominator. 

Solve these for ui and U2 by integration and then plug them 
back into yp to get your particular solution. 

Example 2.4.1. Solve 

y" + y = tan (a:). 

soln: The functions yi = cos(a:) and y2 = sin(x) are funda- 
mental solutions with Wronskian iy(cos(x), sin(x)) = 1. The 
Cramer's rule formulas above become: 



det ( ^ sin(a:) \ / cos(x) 

, Y tan(x) cos(x) J , y — sin(x) tan(a;) 

Ui = , U2 - 



Therefore, 



1 ' " 1 



Ui = -— -, U2 = sm(a:). 



cos(x) ' 

Therefore, using methods from integral calculus, mi = — In | tan(a;)+ 
sec(,T)| + sin(x) and U2 = — cos(a;). Using SAGE, this can be 
check as follows: 

SAGE 



sage : integral (-sin (t) "2/cos (t) ,t) 

-log{sin(t) + l)/2 + log{sin(t) - l)/2 + sin(t) 

sage: integral (cos (t) -sec (t) , t) 

sin(t) - log(tan(t) + sec(t)) 

sage: integral (sin (t) , t) 

-cos (t ) 



As you can see, there are other forms the answer can take. The 
particular solution is 



Up = (— In I tan(a;) + sec(a;)| + sin(a;)) cos(a;) + (— cos(a;)) sin(x). 



The homogeneous (or complementary) part of the solution is 

Uh = ci cos(x) + C2 sin(x), 
so the general solution is 

y =yh + yp = ci cos{x) + C2 sin(a:) 

+(— In I t8Ji{x) + sec(a:)| + sin(x)) cos(x) + (— cos(a:)) sin(a:). 

Using SAGE, this can be carried out as follows: 

SAGE 

sage: SR = SymbolicExpressionRing ( ) 
sage: MS = MatrixSpace ( SR, 2, 2) 

sage: W = MS ( [ [cos (t) , sin (t) ] , [dif f (cos (t) , t) , dif f (sin (t) , t) ] ] : 
sage: W 

[ cos (t) sin (t) ] 

[-sin (t) cos (t) ] 

sage : det (W) 

sin (t) "2 + cos (t) "2 

sage: Ul = MS ( [ [0, sin (t) ] , [tan (t) , dif f (sin (t) , t)]]) 
sage: U2 = MS ( [ [cos (t) , 0] , [diff (cos (t) , t) ,tan (t) ] ] ) 
sage: integral (det (Ul) /det (W) , t ) 

-log(sin(t) + l)/2 + log(sin(t) - l)/2 + sin(t) 
sage: integral (det (U2 ) /det (W) , t ) 
-cos (t) 



Exercise: Use SAGE to solve y" + y = cot(x). 



2.5 Applications of DEs: Spring problems 



2.5.1 Part 1 

Ut tensio, sic vis^. 

- Robert Hooke, 1678 

One of the ways DEs arise is by means of modeling physical 
phenomenon, such as spring equations. For these problems, con- 
sider a spring suspended from a ceiling. We shall consider three 
cases: (1) no mass is attached at the end of the spring, (2) a 
mass is attached and the system is in the rest position, (3) a 
mass is attached and the mass has been displaced fomr the rest 
position. 




Figure 2.5: A spring 
at rest, without mass 
attached. 



Figure 2.6: A spring 
at rest, with mass at- 
tached. 



Figure 2.7: A spring in 
motion. 



■As the extension, so the force. 



One can also align the springs left-to-right instead of top-to- 
bottom, without changing the discussion below. 

Notation: Consider the first two situations above: (a) a spring 
at rest, without mass attached and (b) a spring at rest, with 
mass attached. The distance the mass pulls the spring down is 
sometimes called the "stretch", and denoted s. (A formula for 
s will be given later.) 

Now place the mass in motion by imparting some initial ve- 
locity (tapping it upwards with a hammer, say, and start your 
timer). Consider the second two situations above: (a) a spring 
at rest, with mass attached and (b) a spring in motion. The 
difference between these two positions at time t is called the 
displacement and is denoted x{t). Signs here will be choosen so 
that down is positive. 

Assume exactly three forces act: 

1. the restoring force of the spring, Fgpring, 

2. an external force (driving the ceiling up and down, but may 
be 0), Fexu 

3. a damping force (imagining the spring immersed in oil or 
that it is in fact a shock absorber on a car), F^amp- 

In other words, the total force is given by 

Ftotal Fspring ~l~ Fgxt ~l~ F^Q^^p. 

Physics tells us that the following are approximately true: 

1. (Hooke's law [H-intro]): Fspring = —kx, for some "spring 
constant" A; > 0, 



2. Fext = F{t), for some (possibly zero) function F, 

3. Fdamp = —bv, for some "damping constant" 6 > (where v 
denotes velocity), 

4. (Newton's 2nd law [N-mech]): Ftotai = fnci (where a denotes 
acceleration) . 

Putting this all together, we obtain mx" = ma = —kx + F{t) — 
bv = —kx + F{t) — bx', or 



mx" + bx' + kx = F{t). 

This is the spring equation. When b = F{t) = this is also 
called the equation for simple harmonic motion. 

Consider again first two figures above: (a) a spring at rest, 
without mass attached and (b) a spring at rest, with mass at- 
tached. The mass in the second figure is at rest, so the gravita- 
tional force on the mass, mg, is balanced by the restoring force 
of the spring: mg = ks, where s is the stretch.In particular, the 
spring constant can be computed from the stretch: 




Example: 

A spring at rest is suspended from the ceiling without mass. A 
2 kg weight is then attached to this spring, stretching it 9.8 cm. 
From a position 2/3 m above equilibrium the weight is give a 
downward velocity of 5 m/s. 

(a) Find the equation of motion. 



(b) What is the amplitude and period of motion? 

(c) At what time does the mass first cross equilibrium? 

(d) At what time is the mass first exactly 1/2 m below equilib- 
rium? 

We shall solve this problem using SAGE below. Note m = 2, 
b = F{t) = (since no damping or external force is even 
mentioned), and k = mg/s = 2 • 9.8/(0.098) = 200. There- 
fore, the DE is 2x" + 200x = 0. This has general solution 
x{t) = cicos(lOt) + C2sin(10t). The constants ci and C2 can 
be computed from the initial conditions x{0) = —2/3 (down is 
positive, up is negative), x'{0) = 5. 

Using SAGE , the displacement can be computed as follows: 

SAGE 

sage : t = var ( ' t ' ) 

sage: x = f unction (' x' , t) 

sage: m = var('m') 

sage : b = var ( ' b' ) 

sage: k = var('k') 

sage: F = var('F') 

sage: de = lambda y: m*dif f (y , t , t ) + b*diff(y,t) + k*y - F 
sage : de (x (t ) ) 

-F + m*dif f (x (t) , t, 2) + b*diff(x(t), t, 1) + k*x(t) 
sage: m = 2; b = 0; k = 2*9.8/(0.098); F = 

sage : de (x (t ) ) 

2*diff(x(t), t, 2) + 200. O00000000000*x(t) 
sage : desolve (de (x (t ) ) , [x, t ] ) 
' %kl*sin (10*t) +%k2*cos (10*t) ' 

sage: desolve_laplace (de (x (t) ) , ["t", "x"] , [0,-2/3,5] ) 
' sin (10*t) /2-2*cos (10*t) /3' 



Now we write this in the more compact and useful form A sm{ut-\- 
(j)) using the formulas 



ci cos{ut) + C2 sm{ujt) = A sm{ujt + 0), 
where A = -\/cf~+c| denotes the amplitude and = 2 arctan( ^y|^^ ). 

SAGE 

sage: A = sqrt ( (-2/3) '2+ (1/2) '2) 

sage: A 

5/6 

sage: phi = 2*atan ( (-2/3) / (1/2 + A)) 

sage: phi 

-2*atan (1/2) 

sage: RR(phi) 

-0 . 927295218001612 

sage: sol = lambda t: sin (10*t) /2-2*cos (10*t) /3 
sage: sol2 = lambda t: A*sin(10*t + phi) 
sage: P = plot (sol (t) , 0, 2) 
sage: show(P) 



This is displayed below^. 



^Yoii can also, if you want, type show(plot(sol2(t) ,0,2)) to check that these two 
functions are indeed the same. 



Figure 2.8: Plot of 2x" + 200x = 0, x(0) = 



-2/3, x'{Q) = 5, for < t < 2. 



Of course, the period is 27r/10 = 7r/5 « 0.628. 
To answer (c) and (d), we solve x{t) =0 and x{t) = 1/2: 

SAGE 

sage: solve (A*sin (10*t + phi) == 0,t) 
[t == atan (1/2) /5] 
sage: RR ( atan ( 1 /2 ) / 5 ) 
.0927295218001612 

sage: solve (A*sin (10*t + phi) == 1/2, t) 
[t == (asin{3/5) + 2 *atan ( 1 /2 ) ) /1 ] 
sage: RR ( (asin (3/5) + 2*atan (1/2) ) /lO) 
. 157079632679490 



In other words, a:(0.0927...) « 0, x(0.157...) « 1/2. 



Exercise: Use the problem above. 

(a) At what time does the weight pass through the equilibrium 
position heading down for the 2nd time? 

(b) At what time is the weight exactly 5/12 m below equilibrium 
and heading up? 



2.5.2 Part 2 



Recall from part 1, the spring equation 

mx" + hx -\- kx = F{t) 
where x{t) denotes the displacement at time t. 

Unless otherwise stated, we assume there is no external force: 
Fit) = 0. 

The roots of the characteristic polynomial mD^ + bD = k = 
are 

-b ± - 4m/c 
2m 

There are three cases: 

(a) real distinct roots: in this case the discriminant b^ — Amk 
is positive, so 6^ > 4mA;. In other words, b is "large". This 
case is referred to as overdamped. In this case, the roots 
are negative. 



-b - Vb'^ - ^rnk ^ , -b + - 4mA; 

ri = < 0, and n = < 0, 

2m 2m 

so the solution x{t) = cie^^* + C2e''^* is exponentially de- 
creasing. 

(b) real repeated roots: in this case the discriminant b^ — Amk 
is zero, so 6 = \/Amk. This case is referred to as criti- 
cally damped. This case is said to model new suspension 
systems in cars [D-spr]. 



(c) Complex roots: in this case the discriminant — 4mk is 
negative, so 6^ < Amk. In other words, b is "smah". This 
case is referred to as underdamped (or simple harmonic 
when 6 = 0). 



Example: An 8 lb weight stretches a spring 2 ft. Assume a 
damping force numerically equal to 2 times the instantaneous 
velocity acts. Find the displacement at time t, provided that it 
is released from the equilibrium position with an upward velocity 
of 3 ft/s. Find the equation of motion and classify the behaviour. 

We know m = 8/32 = 1/4, b = 2, k = mg/s = 8/2 = 4, 
x{0) = 0, and x'{0) = —3. This means we must solve 

^x" + 2x' + 4x = 0, x{0) = 0, x'{0) = -3. 

The roots of the characteristic polynomial are —4 and —4 (so 
we are in the repeated real roots case), so the general solution 
is x(t) = cie~^^ + C2te~^*. The initial conditions imply ci = 0, 
C2 = —3, so 

x{t) = -3te-^*. 
Using SAGE , we can compute this as well: 

SAGE 

sage: t = var{'"'"t'') 
sage: x = f unction ''x' ' ) 

sage: de = lambda y: ( 1/4 ) *dif f (y, t , t ) + 2*diff(y,t) + 4*y 

sage : de (x (t ) ) 

diff{x(t), t, 2)/4 + 2*diff(x(t), t, 1) + 4*x{t) 
sage: desolve (de (x (t) ) , [x, t] ) 
' (%k2*t+%kl) *%e'' - (4*t) ' 

sage: desolve_laplace (de (x (t) ) , [ ' 't" , "x' ' ] , [0,0,-3]) 



' -3*t*%e''- (4*t) ' 

sage: f = lambda t : -3*t*e" (-4*t) 
sage: P = plot(f,0,2) 
sage: show(P) 



The graph is shown below. 



Figure 2.9: Plot of (l/4)x"+2x'+4x = 0, x(0) = 0, x'(0) = -3, for < t < 2. 

Example: An 2 kg weight is attached to a spring having spring 
constant 10. Assume a damping force numericaUy equal to 4 
times the instantaneous velocity acts. Find the displacement at 
time provided that it is released from 1 m below equilibrium 
with an upward velocity of 1 ft/s. Find the equation of motion 
and classify the behaviour. 

Using SAGE , we can compute this as well: 



0.75 - 



0.5 - 




-0 75 - 



-] - 



SAGE 



sage : t = var ( ''t ' ' ) 
sage: x = function ( ' 



x' ' ) 



sage: de = lambda y: 2*dif f (y, t, t) + 4*diff(y,t) + 10*y 
sage : desolve_laplace (de (x (t) ) , [ "t", "x"] , [0,1,1]) 
' %e''-t* (sin (2*t) +COS (2*t) ) ' 

sage : desolve_laplace (de (x (t) ) , [ "t", "x"] , [0,1,-1]) 
' %e"-t*cos (2*t) ' 

sage: sol = lambda t: e" (-t) *cos (2*t) 
sage: P = plot ( sol (t ) , , 2 ) 
sage: show(P) 

sage: P = plot ( sol (t ) , , 4 ) 
sage: sliow(P) 



The graph is shown below. 




Figure 2.10: Plot of 2x" + 4x' + lOx = 0, x(0) = 1, x'(0) = -1, for < t < 4. 



Exercise: Use the problem above. Use SAGE to find what time 
the weight passes through the equilibrium position heading down 
for the 2nd time. 

Exercise: An 2 kg weight is attached to a spring having spring 
constant 10. Assume a damping force numerically equal to 4 
times the instantaneous velocity acts. Use SAGE to find the dis- 
placement at time t, provided that it is released from 1 m below 
equilibrium (with no initial velocity). 



2.5.3 Part 3 



If the frequency of the driving force of the spring matches the 
frequency of the homogeneous part Xh{t), in other words if 

x" + LU^x = Fq cos(7t), 

satisfies a; = 7 then we say that the spring-mass system is in 
(pure, mechanical) resonance. 

Example: Solve 

x" + u'^x = Fo cos(7t), x{0) = 0, x'{0) = 0, 

where u; = 7 = 2 (ie, mechanical resonance). We use SAGE for 
this: 

SAGE 

sage : t = var ( ' t' ) 

sage: x = f unction (' x' , t) 

sage: (m, b, k, w, FO ) = var ( "m, b, k, w, FO " ) 

sage: de = lambda y: diff(Y,t,t) + vi~2*y - FO*cos(w*t) 



sage: m=l; b=0; k=4; F0=1; w=2 
sage : de solve (de (x (t) ) , [x, t] ) 

' (2*t*sin (2*t) +COS (2*t) ) /8+%kl*sin (2*t) +%k2*cos (2*t) ' 
sage: desolve_laplace (de (x (t ) ) , ["t", "x"] , [0,0,0] ) 

't*sin (2*t) /4' 
sage: soln = lambda t : t*sin(2*t)/4 
sage: P = plot (soln (t) , 0, 10) 
sage: show(P) 



This is displayed below: 




Figure 2.11: A forced undamped spring, with resonance. 



Example: Solve 

+ u^x = Fo cos(7t), x{0) = 0, x'{0) = 0, 

where cj = 2 and 7 = 3 (ie, mechanical resonance). We use 
SAGE for this: 

SAGE 

sage : t = var ( ' t ' ) 

sage: x = function (' x' , t) 

sage: (m, b, k, w, g, FO ) = var ( "m, b, k, w, g, FO " ) 
sage: de = lambda y: diff(y,t,t) + w"2*y - FO*cos(g*t) 
sage: m=l; b=0; k=4; F0=1; w=2; g=3 
sage: desolve_laplace (de (x (t) ) , ["t", "x" ] , [0,0,0]) 

' cos (2*t) /5-cos (3*t) /5' 
sage: soln = lambda t : cos (2*t) /5-cos (3*t) /5 
sage: P = plot (soln (t) , 0, 10) 
sage: show(P) 



This is displayed below: 



Figure 2.12: A forced undamped spring, no resonance. 



2.6 Applications to simple LRC circuits 



An LRC circuit is a closed loop containing an inductor of L hen- 
ries, a resistor of R ohms, a capacitor of C farads, and an EMF 
(electro-motive force), or battery, of E{t) volts, all connected in 
series. 

They arise in several engineering applications. For example, 
AM/FM radios with analog tuners typically use an LRC circuit 
to tune a radio frequency. Most commonly a variable capaci- 
tor is attached to the tuning knob, which allows you to change 
the value of C in the circuit and tune to stations on different 
frequencies [R-cir]. 

We use the following "dictionary" to translate between the 
diagram and the DEs. 



bjbj object 


term m Dh, 


units 


■u 1 

symbol 


charge 


q = f i{t) dt 


coulombs 




current 


i = q' 


amps 




emf 


e = e{t) 


volts V 




resistor 


Rq' = Ri 


ohms Q 




capacitor 




farads 


— 1 1— 


inductor 


Lq" = Li' 


henries 





Kirchojf's First Law: The algebraic sum of the currents trav- 
elling into any node is zero. 

Kirchoff 's Second Law: The algebraic sum of the voltage drops 
around any closed loop is zero. 

Generally, the charge at time t on the capacitor, satisfies 
the DE 

Lq" + Rq' + ^q = Eit). (2.7) 

Example 1: Consider the simple LC circuit given by the fol- 
lowing diagram. 

1 henry 
sin(2t)+sin(11t) volts 

< 1 I— < ^ 

1/C farads 

Figure 2.13: A simple LC circuit. 



According to Kirchoff's 2"^ Law and the above "dictionary", 
this circuit corresponds to the DE 

g" + ig = sin(2t) + sin(llt). 
The homogeneous part of the solution is 



qh{t) = ci cos{t/VC) + ci sin(t/VC). 
If C ^ 1/4 and C ^ 1/121 then 

Qpit) = ^}^_ 4 sin(2t) + \ sin(llt). 

When C = 1/4 and the initial charge and current are both 
zero, the solution is 

qii) = -^sin(llt) + i^sin(2i) - iicos(2i). 

, SAGE 



sage : t = var ( "t " ) 

sage: q = function ( "q" , t) 

sage: L,R,C = var("L,R,C") 

sage: E = lambda t : sin (2*t) +sin (ll*t) 

sage: de = lambda y: L*dif f (y, t , t ) + R*diff(y,t) + (1/C)*y-E(t) 

sage: L, R, C=l, 0, 1/4 

sage : de (q (t ) ) 

diff{q(t), t, 2) - sin(ll*t) - sin{2*t) + 4*q(t) 

sage: desolve_laplace (de {q(t) ) , ["t", "q"] , [0, 0,0] ) 
' -sin (ll*t) /117+161*sin (2*t) /93 6-t*cos (2*t) /4' 

sage: soln = lambda t: -sin (ll*t) /117+161*sin {2*t) /936-t*cos (2*t 

sage: P = plot (soln, 0, 10) 

sage: show(P) 



/4 



This is displayed below: 



Figure 2.14: A LC circuit, with resonance. 



You can see how the frequency cj = 2 dominates the other 
terms. 

When < R < LjC the homogeneous form of the charge 
in (2.7) has the form 

q^{t) = cie"*cos(/?t) + C2e"^ sin(/3t), 

where a = -R/2L < and = ^JAL/C - / {2L). This is 
sometimes caUed the transient part of the solution. The re- 
maining terms in the charge are called the steady state terms. 

Example: An LRC circuit has a 1 henry inductor, a 2 ohm 
resistor, 1/5 farad capacitor, and an EMF of 50cos(t). If the 
initial charge and current is 0, since the charge at time t. 
The I VP describing the charge q{t) is 

q" + 2q + 5g = 50cos(i), g(0) = g'(0) = 0. 
The homogeneous part of the solution is 



Qhit) = cie *cos(2t) +C2e *sin(2t). 

The general form of the particular solution using the method of 
undetermined coefficients is 



sage: t = var("t") 

sage: q = f unction ( "q" , t ) 

sage: L,R,C = var("L,R,C") 

sage: E = lambda t: 50*cos(t) 

sage: de = lambda y: L*dif f (y, t, t) + R*diff(y,t) + (1/C)*y-E(t) 

sage: L,R,C = 1,2,1/5 

sage : de (q (t ) ) 

diff(q(t), t, 2) + 2*diff(q(t), t, 1) + 5*q(t) - 50*cos(t) 

sage: desolve_laplace (de (q(t) ) , ["t", "q"] , [0,0,0]) 
' %e''-t* (-15*sin (2*t) /2-10*cos (2*t) ) +5*sin (t) +10*cos (t) ' 

sage: soln = lambda t:\ 

e'~ (-t) * (-15*sin (2*t) /2-10*cos (2*t) ) +5* sin (t) +10*cos (t) 

sage: P = plot (soln, 0, 10) 

sage: show(P) 

sage: P = plot (soln, 0, 20) 

sage: show(P) 

sage: soln_ss = lambda t: 5*sin (t ) +10*cos (t ) 

sage: P_ss = plot (soln_ss, 0, 10) 

sage: soln_tr = lambda t: e' (-t) * (-15*sin (2*t) /2-10*cos (2*t) ) 

sage: P_tr = plot (soln_tr, 0, 10, linestyle=" — ") 

sage: show(P+P_tr) 



This plot (the solution superimposed with the transient part of 
the solution) is displayed below: 



qp{t) = Ai cos(t) + A2 sin(t). 



Solving for Ai and A2 gives 




SAGE 



Figure 2.15: A LRC circuit, with damping, and the transient part (dashed) of the 
solution. 

Exercise: Use SAGE to solve 

q +^q = sm{2t) + sin(lU), q{0) = q\0) = 0, 
in tiie case C = 1/121. 



2.7 The power series method 
2.7.1 Part 1 



In this part, we recall some basic facts about power series and 
Taylor series. We will turn to solving DEs in part II. 

Roughly speaking, power series are simply infinite degree poly- 
nomials 

oo 

f{x) = flo + CLix + a2X^ + ... = ttkx'^, (2.8) 

k=0 

for some real or complex numbers ao,ai,... The number is 
called the coefficient of x^, for k = 0, 1, .... Let us ignore for 
the moment the precise meaning of this infinite sum (How do 
you associate a value to an infinite sum? Does the sum converge 
for some values of xl If so, for which values? ...) We will return 
to that later. 

First, some motivation. Why study these? This type of func- 
tion is convenient for several reasons 

• it is easy to differentiate (term-by-term): 

oo oo 

f'{x) = ai-\-2a2X-\-3a3x'^-\-... = ^^/cafcx'^""^ = ^^(A;+l)afc+ix^, 

A;=0 k=0 

• it is easy to integrate (term- by-term): 



f 11 oo ^ oo ^ 

/ f{x) dx = aQX-\--aix'^-\--a2X^-\-... = ^ ^rqr^;"^;^^^^^ = ^ -j^CLk+iX 

A;=0 k=l 



if (as is often the case) the a^s tend to zero very quickly, 
then the sum of the first few terms of the series are often a 
good numerical approximation for the function itself, 

power series enable one to reduce the solution of certain dif- 
ferential equations down to (often the much easier problem 
of) solving certain recurrance relations. 

Power series expansions arise naturally in Taylor's Theorem 
of the Mean: If f{x) is n + 1 times continuously differen- 
tiable in {a,x) then there exists a point ^ G {a,x) such 
that 



fix) = f(a) + (x- a}f(a) + ^^^/"(a) + • • • 

+ ^-^^f^"\a) + <-;7f"' /<"""(4). (2.9) 
n! (n + Ij! 

The sum 



Tnix) = f{a)+(x-a)na) + ^-^^f"{a)+- ■ 

is called the n-th degree Taylor polynomial of / cen- 
tered at a. For the case n = 0, the formula is 

fix) = fia) + ix-a)f'iO, 

which is just a rearrangement of the terms in the theorem 
of the mean, 

X — a 



Some examples: 



• Geometric series: 



1 — X 



oo 



= ^x" (2.10) 



n=0 



To see this, assume \x\ < 1 and let n ^ oo in the polyno- 
mial identity 

1 + X + H h x"""^ = . 

1 — X 

For X > 1, the series does not converge. 



• The exponential function: 



2 3 4 

2 3 4 

tXy tAj tXy 



n=0 



To see this, take f{x) = e'^ and a = in Taylor's theorem 



(2.9), using the fact that 4^e^' = and = 1: 



X2 /y*3 /v>^ 

e" = l + a; + — + — + ••• + — + 



2! 3! n\ (n + 1)! 



for some ^ between and x. Perhaps it is not clear to 
everyone that as n becomes larger and larger {x fixed) , the 
last ("remainder") term in this sum goes to 0. However, 
Stirling's formula tells us how large the factorial function 
grows, 



n!~v^ - (l + 0(-)), 
Ve/ n 

so we may indeed take the limit as n ^ oo to get (2.11). 

Wikipedia's entry on "Power series" [Pl-ps] has a nice an- 
imation showing how more and more terms in the Taylor 
polynomials approximate better and better. 



• The cosine function: 



2 4 fi 

rp^ rp^ rf^ 



COS X = \ — ^ -\- iri — + 
2 24 720 

2 4 fi 

nf*^ rp^ 

Ju tXj tAJ 

~2[^4[~6[^" 



= E(-i)"(^ (2-12) 

n=0 ^ ^ 



This too follows from Taylor's theorem (take f{x) = cosx 
and a = 0). However, there is another trick: Replace x in 
(2.11) by ix and use the fact ("Euler's formula") that e*^ = 
cos(x) + zsin(x). Taking real parts gives (2.12). Taking 
imaginary parts gives (2.13), below. 



• The sine function: 



^ 7 



sm X = X 1 h • • • 

6 120 5040 

^ 7 

tty tAJ Ju 

~ 3! 5! ~ tT ^ " " " 

00 ™2n+l 
= y(-ir— (2.13) 

^ 2n+l! ^ ^ 

n=0 ^ ^ 

Indeed, you can formally check (using formal term-by-term 
differentiation) that 

— -^cos(a;) = sin(a:). 
ax 

(Alternatively, you can use this fact to deduce (2.13) from 
(2.12).) 

• The logarithm function: 

log(l — x) = —X x^ x^ — -x"^ + • • • 

&\ J 2 3 4 

OO 

= -E„^" (2-14) 

n=0 

This follows from (2.10) since (using formal term-by-term 
integration) 

rx ^ 

SAGE 



sage: taylor (sin (x) , x, 0, 5) 
X - x"3/6 + x~5/120 



sage: PI = plot ( sin (x) , , pi ) 

sage: P2 = plot (x, , pi , linestyle=" — ") 

sage: P3 = plot (x-x" 3/ 6 , , pi, linestyle=" - . " ) 

sage: P4 = plot (x-x" 3/ 6+x" 5 /12 , , pi, linestyle=" : " ) 

sage: Tl = textC'x", (3,2.5)) 

sage: T2 = text ( "x-x" 3/3 ! " , ( 3 . 5 , -1) ) 

sage: T3 = text ( "x-x" 3/3 ! +x" 5/ 5 ! " , ( 3 . 7 , . 8 ) ) 

sage: T4 = text ( " sin (x) " , ( 3 . 4 , . 1 ) ) 

sage: show (P1+P2+P3+P4+T1+T2+T3+T4 ) 



This is displayed below: 




Figure 2.16: Taylor polynomial approximations for sin(x) 



Exercise: Use SAGE to plot successive Taylor polynomial ap- 
proximations for cos(a:). 

Finally, we turn to the meaning of these sums. How do you 
associate a value to an infinite sum? Does the sum converge for 
some values of xl If so, for which values? . We will (for the 
most part) answer all of these. 

First, consider our infinite power series f{x) in (2.8), where the 
ttk are all given and x is fixed for the momemnt. The partial 



sums of this series are 



/o(x) = ao, fi{x) = ao + aix, /2(x) = ao + aix + a2X^, • • • • 

We say that the series in (2.8) converges at x if the hmit of 
partial sums 

hm fn{x) 

exists. There are several tests for determining whether or not a 
series converges. One of the most commonly used tests is the 

Root test: Assume 

L = lim \akx''\'^/^ = \x\ lim \ak\^^'' 

exists. If L < 1 then the infinite power series f{x) in (2.8) 
converges at x. In general, (2.8) converges for all x satisfying 

- lim \ak\~^^'' < X < lim \ak\'^^''. 

k^oo k^oo 

The number limfc^oo |flfc| ^'^^ (if it exists, though it can be oo) is 
called the radius of convergence. 

Example: The radius of convergence of (and cos(a:) and 
sin(x)) is oo. The radius of convergence of 1/(1 — x) (and log(H- 
x)) is 1. 

Example: The radius of convergence of 

fc=0 

can be determined with the help of SAGE . We want to compute 



SAGE 



sage : k = var { ' k' ) 

sage: limit ( ( (k'T+k+l) / (2'k+k'2) ) ' (-1/k) ,k=infinity) 
2 



In other words, the series converges for all x satisfying — 2 < 
X <2. 

Exercise: Use SAGE to find the radius of convergence of 

/(-) = E 

A;=0 



2.7.2 Part 2 



In this part, we solve some DEs using power series. 
We want to solve a problem of the form 

y"{x) + p{x)y'{x) + y{x) = /(x), (2.15) 

in the case where p{x), q{x) and f{x) have a power series expan- 
sion. We will call a power series solution a series expansion 
for y{x) where we have produced some algorithm or rule which 
enables us to compute as many of its coefficients as we like. 



Solution strategy: Write y{x) = ao + clix + a2x'^ + ... = 
Ylh=o^kx'^^ some real or complex numbers ao,ai, .... 

• Plug the power series expansions for y, p, q, and / into the 
DE (2.15). 

• Comparing coeffiients of like powers of x, derive relations 
between the afs. 

• Using these recurrance relations [R-ps] and the ICs, solve 
for the coefficients of the power series of y{x). 

Example: Solve y' — y = 5, y{0) = —4, using the power series 
method. 

This is easy to solve by undetermined coefficients: yh{x) = cie^ 
and yp{x) = Ai. Solving for Ai gives Ai = —5 and then solving 
for ci gives —4 = y{0) = — 5 + cie° so ci = 1 so y = — 5. 

Solving this using power series, we compute 

y'{x) = ai + 2a2a^ + Sasx^ + ... = X]^q(A; + l)aA;+ia;^ 
-y{x) = -ao - aix - a2X^ - ... = J2T=o~^kx'' 

5 = (-ao + ai) + (-ai + 2a2)x + ... = Xl^o(~^fc + + 

Comparing coefficients, 

• for A; = 0: 5 = — ao + ai, 

• for A; = 1: = — ai + 2a2, 

• for general k: = —ak + (A; + l)afc+i for A; > 0. 



The IC gives us —4 = y(0) = ao, so 

ao = -4, tti = 1, a2 = 1/2, aa = 1/6, ••• ,afc = l/A;!. 
This implies 

y{x) = + X + +x^/k\+ =-5 + e^, 

which is in agreement from the previous discussion. 

Example: Solve Bessel's equation [B-ps] of the 0-th order, 

X V + ^y' + = 0, ?/(0) = l, ?/'(0) = 0, 

using the power series method. 

This DE is so well-known (it has important applications to 
physics and engineering) that the series expansion has already 
been worked out (most texts on special functions or differential 
equations have this but an online reference is [B-ps]). Its Taylor 
series expansion around is: 

m=0 

for all X. We shall see below that y{x) = Jo(a:). 

Let us try solving this ourselves using the power series method. 
We compute 



x'^y"{x) 

xy'{x) 

x^y{x) 





= + • X + 2a2x'^ + Gasx^ + 12a4X^ + ... = EfcLo(^ + 2)(A; + 
= + aix + 2a2x'^ + Sa-^x^ + ... = YlT=o kakX^ 
= + ■ x + aox^ + aix^ + ... = Yl'k=2 (^k-2X^ 

. = aix + Yjk=2i^k-2 + k'^ak)x^. 



= + aix + (ao + 4:a2)x^ + 



By the ICs, ao = 1, ai = 0. Comparing coefficients, 



k'^ttk = -afc-2, A; > 2, 

which imphes 

1 11 111 

a2 = -{^f, 03 = 0, a4 = (-•-)^ a5 = 0, as = -(-•-•-)^ • • • . 

In general, 

a2k = (-1)^2-2^-^^, a2k+i = 0, 

for /c > 1. This is in agreement with the series above for Jq. 

Some of this computation can be formally done in SAG Busing 
power series rings. 

SAGE 

sage: R6 . <aO , al , a2 , a3 , a4 , a5 , a6> = PolynomialRing (QQ, 7 ) 
sage: R.<x> = PowerSeriesRing (R6) 

sage: y = aO + al*x + a2*x'~2 + a3*x"3 + a4*x'4 + a5*x'5 +\ 

a6*x'6 + 0(x"7) 
sage: yl = y . derivative ( ) 
sage: y2 = yl . derivative ( ) 

sage: x"2*y2 + x*yl + x'"2*y 

al*x + (aO + 4*a2)*x"2 + (al + 9*a3)*x"3 + (a2 + 16*a4)*x^4 +\ 
(a3 + 25*a5)*x"5 + (a4 + 36*a6)*x"6 + 0(x"7) 



This is consistent with our "paper and pencil" computations 
above. 

SAGE knows quite a few special functions, such as the various 
types of Bessel functions. 

SAGE 

sage: b = lambda x:bessel_J (x, 0) 



sage: P = plot (b, 0, 20, thickness=l) 

sage: show(P) 

sage: y = lambda x: 1 - (l/2)'~2*x"2 + (l/8)~2*x"4 - (l/48)"2*x"6 

sage: PI = plot (y, 0, 4, thickness=l) 

sage: P2 = plot (b, , 4 , linestyle=" — ") 

sage: show(Pl+P2) 



This is displayed below: 




Figure 2.17: The Bessel function Figure 2.18: A Taylor polynomial 
Jo{x), for < a; < 20. approximation for Jo{x). 



Exercise: Use SAGE to find the first 5 terms in the power series 
solution to y" -\-y = 0, y{0) = 1, y'{0) = 0. Plot this Taylor 
polynomial approximation over — tt < x < it. 



2.8 The Laplace transform method 



2.8.1 Part 1 

The Laplace transform (LT) of a function f{t), defined for all 
real numbers t > 0, is the function -F(s), defined by: 



This is named for Pierre-Simon Laplace, one of the best French 
mathematicians in the mid-to-late 18th century [L-lt], [LT-lt]. 
The LT sends "nice" functions of t (we will be more precise later) 
to functions of another variable s. It has the wonderful property 
that it transforms constant-coefficient differential equations in t 
to algebraic questions in s. 

The LT has two very familiar properties: Just as the integral 
of a sum is the sum of the integrals, the Laplace transform of a 
sum is the sum of Laplace transforms: 



Just as constant factor can be taken outside of an integral, the 
LT of a constant times a function is that constant times the LT 
of that function: 



In other words, the LT is linear. 

For which functions / is the LT actually defined on? We want 
the indefinite integral to converge, of course. A function f{t) is 
of exponential order a if there exist constants to and M such 




C [fit] + 9{t)] = C [fit]] + C [g{t)] 



C [afit]] = aC [fit]] 



that 



If Jq° f{t) dt exists and f{t) is of exponential order a then the 
Laplace transform C [/] (s) exists for s > a. 



\f{t)\ <Me"*, for alH>to- 



Example 2.8.1. Consider the Laplace transform of f{t) 
The LT integral converges for s > 0. 



= 1. 



/•OO 

jC[f]{s)= / e-^'dt 
Jo 



— e 
s 



-St 



OO 



1 

S 



Example 2.8.2. Consider the Laplace transform of fit) 
The LT integral converges for s > a. 



poo 

£[/](s)= / e^^'-'^^dt 
Jo 



s — a 



,{a-s)t 



s — a 



Example 2.8.3. Consider the Laplace transform of the unit step 
(Heaviside) function, 



u{t — c) = 



for t < c 

1 for t > c, 



where c> 0. 



POO 

C[u{t -c)]= / e-'*H{t - c) dt 
Jo 



'0 

POO 

L " 



c 



— S 
-cs 



oo 



= for s > 

s 

The inverse Laplace transform in denoted 

fit) = c-'[Fism), 

where F{s) = C[fit)] (s). 
Example 2.8.4. Consider 



1, Vt<2, 
0, ont> 2. 

We show how SAGE can be used to compute the LT of this. 

SAGE 

sage : t = var ( ' t' ) 

sage : s = var ( ' s' ) 

sage: f = Piecewise ([[ (0, 2) , 1] ,[ (2, infinity) , 0] ] ) 

sage: f.laplace(t, s) 

1/s - e' (- (2*s) ) /s 

sage: fl = lambda t: 1 

sage: f2 = lambda t: 

sage: f = Piecewise ([[ (0, 2) , fl] ,[ (2, 10) , f2] ] ) 

sage: P = f .plot (rgbcolor= (0 . 7, . 1, . 5) , thickness=3) 

sage: show(P) 



According to SAGE, jC [f] (s) = 1/s — e~^*/s. Note the function 
f was redefined for plotting purposes only (the fact that it was 
redefined over < ^ < 10 means that SAGE will plot it over that 
range.) The plot of this function is displayed below: 



Figure 2.19: The piecewise constant function 1 — u{t — 2). 

Next, some properties of the LT. 
• Differentiate tfie definition of the LT with respect to s: 

■'OO 



POO 

F'{s) = - e-'Hf{t) dt. 
Jo 



Repeating this: 



"OO 



;F{s) = {-iri e-^H-f{t)dt 



(2.16) 



• In the definition of the LT, replace f{t) by it's derivative 
/'(«): 



poo 

^[nt)]{s)= / e-^'f{t)dt. 
Jo 



Now integrate by parts {u = e dv = fit) dt): 

poo noo 

/ e-^'f{t)dt = f{t)e-^X- f{t)i-sye-^'dt = -fiO)+sC[f{t)] 
Jo Jo 

Therefore, if F{s) is the LT of f{t) then sF{s) - /(O) is the 
LTof fit): 

c[nt)]is) = sc[m]is)-m. (2.17) 

• Replace / by f in (2.17), 

c[nt)]is) = sc[nt)]is)-noi (2.18) 

and apply (2.17) again: 

C [/"(*)] (^) = s'C [/Wl (s) - ./(O) - /'(O), (2.19) 

• Using (2.17) and (2.19), the LT of any constant coefficient 
ODE 

ax"{t) + bx'{t) + cx{t) = f{t) 

is 

a{s'^/:[x{t)] {s)-sx{0)-x'{0)) + b{sC [x{t)] (s)-x(0)) + c£ [x{t)] {s) = F{s), 

where F{s) = C[f{t)] (s). In particular, the LT of the 
solution, X{s) = C [x{t)] (s), satisfies 



X{s) = (F(s) + asx{0) + ax'{0) + bx{0))/{as'^ + 6s + c). 



Note that the denominator is the characteristic polynomial 
of the DE. 



Moral of the story: it is always very easy to compute the LT 
of the solution to any constant coefficient non-homogeneous 
linear ODE. 

Example 2.8.5. We know now how to compute not only the 
LT of f{t) = (it's F{s) = (s - a)-^) but also the LT of any 
function of the form t^e^^ by differentiating it: 



C [te«^] = -F'{s) = (s-a)-2, C [t\^'] = F"(s) = 2.(s-a)-^ C [t^e""'] = -F'{s) 
and in general 

C [re"*] = -F'(s) = n! • (s - a)-^-\ (2.20) 
Example 2.8.6. Let us solve the DE 

x' + x = t^^^e-\ x{0) = 0. 

using LTs. Note this would be highly impractical to solve using 
undetermined coefficients. (You would have 101 undetermined 
coefficients to solve for!) 

First, we compute the LT of the solution to the DE. The LT of 
the LHS: by (2.20), 

C[x' + x\ = sX(s)+X(s), 
where F{s) = C [f{t)] (s). For the LT of the RHS, let 

1 



By (2.16), 



jlOO jlOO 

F{s) = C [t'^^'e-'] = ^ 



1 



The first several derivatives of ^ are as follows 



d 1 



1 



d^ 



1 



= 2- 



1 



1 



dss + l (s + l)2' (/s2s + l (s + l)3' ^^s3s + l 
and so on. Therefore, the LT of the RHS is: 

^/loo 1 . 1 



= -62 



is + 1)^ 



^^sloo 8 + 1 



= 100! 



Consequently, 



X{s) = 100!- 



(5 + 1) 



1 



101 



(S + 1)102 • 

Using (2.20), we can compute the ILT of this: 



x{t) = [X(s)] = ^100!^^ ^ 
Example 2.8.7. Let us solve the DE 



1 



101 



101!- 



s + 1) 



102 



101 



+ 2x +2x = a;(0) = a;'(0) = 0, 

M5m^ LTs. 

The LT of the LHS: by (2.20) and (2.18), 

C [x" + 2x' + 2x\ = (s^ + 2s + 2)X(s), 
as in the previous example. The LT of the RHS is: 

1 



C [e-T = 



S + 2" 



Solving for the LT of the solution algebraically: 



^^'^ (s + 2)((s + 1)2 + 1) • 

The inverse LT of this can be obtained from LT tables after 
rewriting this using partial fractions: 



Xis) = 



1 1 



1 



1 1 



1 s + 1 



2 s + 2 2 (s + 1)2 + 1 
The inverse LT is: 



2 s + 2 2 (s + 1)2 + 1 2 (s + 1)2 + 1' 



xit) = [Xis)] = I ■ - 1 • e-*cos(t) + ^ • e-'sm{t). 
We show how SAGE can be used to do some of this. 

SAGE 



sage : t = var { ' t ' ) 

sage : s = var { ' s ' ) 

sage: f = 1/ ( (s+2) * ( (s+1) "2+1) ) 

sage: f .partial_f raction ( ) 

l/(2*(s + 2)) - s/(2*(s'^2 + 2*s + 2)) 
sage: f . inverse_laplace { s , t ) 

e" (-t) * (sin (t) /2 - cos(t)/2) + e"(-(2*t))/2 



Exercise: Use SAGE to solve the DE 



x" + 2x' + bx = e-\ x{0) = x'{0) = 0. 



2.8.2 Part 2 



In this lecture, we shall focus on two other aspects of Laplace 
transforms: 

• solving differential equations involving unit step (Heavi- 
side) functions, 

• convolutions and applications. 

It follows from the definition of the LT that if 

m^Fis) = c[fms), 

then 

f^t)uit - c) ^ e-''C[f{t + c)](s), (2.21) 

and 

f{t - c)u{t - c) ^ e-"'F(s). (2.22) 
These two properties are called translation theorems. 

Example 2.8.8. First, consider the Laplace transform of the 
piecewise- defined function fit) = {t — l)^u{t — 1). Using (2.22), 
this is 

Cm] = e-^C[t\s) = 2\e-\ 

Second, consider the Laplace transform of the piecewise- constant 
function 

!0 fort< 0, 
-1 for0<t<2, 
1 fort> 2. 



This can be expressed as f{t) = 



-u{t) + 2u{t - 2), so 



mt)] 



= -jC[u{t)] + 2jC[u{t-2)] 




Finally, consider the Laplace transform of f{t) = sin(^)w(^— tt). 
Using (2.21), this is 



C[f{t)] = e^'''jC[sm{t+7T)]{s) = e-""' jC[- sm{t)]{s) = -e" 



The plot of this function f{t) = sm{t)u{t — tt) is displayed below: 



Figure 2.20: The piecewise continuous function u{t — vr) sin(t). 
We show how SAGE can be used to compute these LTs. 



sage : t = var ( ' t ' ) 

sage : s = var ( ' s ' ) 

sage: assume (s>0) 

sage: f = Piecewise ([[(0,1), 0], [(1, infinity) , (t-1) "2] ] ) 

sage: f.laplace(t, s) 



1 




10 



SAGE 



2*e" (-S) /s^3 

sage: f = Piecewise ( [ [ (0, 2) , -1] , [ (2, infinity) , 2] ] ) 
sage: f.laplace(t, s) 

3*e" (- (2*s) ) /s - 1/s 

sage : f = Piecewise ( [ [ {0,pi) ,0] , [ (pi, infinity) ,sin(t) ] ] ) 
sage: f.laplace(t, s) 
-e" (- (pi*s) ) / (3-2 + 1) 
sage: fl = lambda t: 
sage: f2 = lambda t: sin(t) 

sage: f = Piecewise ([[( , pi ), fl ],[ (pi, 10 ), f 2 ]] ) 
sage: P = f . plot (rgbcolor= ( . 7 , . 1 , . 5) , thickness=3 ) 
sage: show(P) 



The plot given by these last few commands is displayed above. 

Before turning to differential equations, let us introduce con- 
volutions. 

Let f{t) and g{t) be continuous (for t > - for t < 0, we 
assume f{t) = g{t) = 0). The convolution of f{t) and g{t) is 
defined by 



The LT of the convolution is the product of the LTs. (Or, equiv- 
alently, the inverse LT of the product is the convolution of the 
inverse LTs.) 

To show this, do a change- of- variables in the following double 
integral: 




The convolution theorem states 



Cif * gms) = F{s)G(s) 



Cif](3)Clg]{s). 



oo 

St 



JO Jo 

"OO />oo 



/ 
I 

Jo 



f{u)g{t — u) dudt 

poo poo 

/ / e-'^f{u)g{t-u)dtdu 

Jo Ju 

poo 

e"'VH / e''^*-''^g{t-u)dtdu 

J u 

poo 

e-'^f(u)du / e-"'g(v)dv 
= C[f\(s)C[g\[s). 
Example 2.8.9. Consider the inverse Laplace transform of 



This can he computed using partial fractions and LT tables. 
However, it can also be computed using convolutions. 
First we factor the denominator, as follows 



1 



1 1 



— s — 1 



We know the inverse Laplace transforms of each term: 



1 



1 



s-1 



We apply the convolution theorem: 



1 1 
s — 1 



= / ue* ^du 
Jo 

= e' [-ue-']l-e' f 

JO 

= -t - 1 + e* 



-e du 



Therefore, 



" 1 1 ■ 
s — I 



{t) = e^ -t- 1. 



Example 2.8.10. Here is a neat application of the convolution 
theorem. Consider the convolution 

f{t) = 1*1*1*1*1. 

What is it? Firsts take the LT. Since the LT of the convolution 
is the product of the LTs: 

£[1 * 1 * 1 * 1 * l](s) = (l/s)5 = ^ = F{s). 
We know from LT tables that [^] {t) =t^, so 



m = [F(s)] (t) = 



4! 



Now let us turn to solving a DE of the form 



ay" + hy' + cy = f[t), y{0) = yo, y'{0)=yi. 
First, take LTs of both sides: 



(2.23) 



as^y(s) — asyo — ayi + bsY{s) — byo + cY{s) = F{s). 



so 



Yis) = 



The function 



1 



as^ + 6s + c 



F{s) + 



asyo + ayi + byp 
as'^ + 6s + c 



(2.24) 



as-^+bs+c sometimes called the transfer function 
(this is an engineering term) and it's inverse LT, 



w{t) = 



1 



it), 



as^ + 5s + c 
the weight function for the DE. 

Lemma 2.8.1. If a then w{t) = 0. 

(The only proof I have of this is a case-by-case proof using LT 
tables. Case 1 is when the roots of as^ + 6s + c = are real and 
distinct, case 2 is when the roots are real and repeated, and case 
3 is when the roots are complex. In each case, w{0) = 0. The 
verification of this is left to the reader, if he or she is interested.) 

By the above lemma and the first derivative theorem. 



w'{t) = 



it). 



as^ + 6s + c 

Using this and the convolution theorem, the inverse LT of 
(2.24) is 



y{t) = {w * f ){t) + avQ ■ w'{t) + {ayi + hy^) ■ w{t). (2.25) 
This proves the following fact. 

Theorem 2.8.1. The unique solution to the DE (2.23) is (2.25). 



Example 2.8.11. Consider the DEy"+y = I, y{0) = y'{0) = 1. 

The weight function is the inverse Laplace transform of ^ 
so w{t) = sm{t). By (2.25), 



S2 + 1' 



y[t) = 1 * sin(i) = / sin(w) du = — cos(w)|o = 1 — cos(i). 
(Yes, it is just that easy!) 



If the "impulse" f{t) is piecewise-defined, sometimes the con- 
volution term in the formula (2.25) is awkward to compute. 

Example 2.8.12. Consider the DE y" - y' = u{t - I), y{0) = 
y'iO) = 0. 

Taking Laplace transforms gives s'^Y(s) — sY(s) = -e"*, so 



Yis) = 



1 



— s-^ 



We know from a previous example that 

1 



— 



(t) = e*-t-l, 



so by the translation theorem (2.22), we have 



yit) = 



— 



(t) = (e*-i-(t-l)-l).w(t-l) = (e'-'-tyu{t-l) 



t-i 



At this stage, SAGE lacks the functionality to solve this DE. 



Exercise: (a) Use SAGE to take the LT of u{t — 7r/4) cos(t). 
(b) Use SAGE to compute the convolution sin(t) * cos(t). 



Chapter 3 

Systems of first order 
differential equations 
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3.1 An introduction to systems of DEs: Lanch- 
ester's equations 



The goal of military analysis is a means of reliably predicting 
the outcome of military encounters, given some basic informa- 
tion about the forces' status. The case of two combatants in an 
"aimed fire" battle was solved during World War I by Frederick 
William Lanchester, a British engineer in the Royal Air Force, 
who discovered a way to model battle-field casualties using sys- 
tems of differential equations. He assumed that if two armies 
fight, with x{t) troops on one side and y{t) on the other, the 
rate at which soldiers in one army are put out of action is pro- 
portional to the troop strength of their enemy. This give rise to 
the system of differential equations 

r x\t) = -Ay{t), xiO)=xo, 
\ y'{t) = -Bx{t), 2/(0) = yo, 

where A > and 5 > are constants (called their fighting effec- 
tiveness coefficients) and xq and are the intial troop strengths. 
For some historical examples of actual battles modeled using 
Lanchester 's equations, please see references in the paper by 
McKay [M-intro]. 

We show here how to solve these using Laplace transforms. 

Example: A battle is modeled by 

( x' = -4y, x{0) = 150, 
\y' = -X, yiO) = 90. 

(a) Write the solutions in parameteric form, (b) Who wins? 
When? State the losses for each side. 



soln: Take Laplace transforms of both sides: 



sL[x it)] (s) -a;(0) 



4L[y(t)] (s) 



sL [x {t)] (s) - X (0) 
Solving these equations gives 



AL[y{t)] {s). 



L\yit)] is) 



L [x (t)] is) 



sx (0) -4y(0) _ 150s -360 
s2-4 ~ s2-4 ' 

-sy (0) + X (0) _ -90 s + 150 
s2 - 4 s2 - 4 



Laplace transform Tables give 

x{t) = -15e^* + 165e-2* 

yit) = 90 cosh (2 1) - 75 sinh (2 
Their graph looks like 

The "?/-army" wins. Solving for x{t) = gives = log(ll)/4 = 
.5994738182..., so the number of survivors is y{t^^n) = 49.7493718, 
so 49 survive. 

Lanchester's square law. Suppose that if you are more inter- 
ested in y as a function of x, instead of x and y as functions of 
t. One can use the chain rule form calculus to derive from the 
system x'{t) = —Ay{t),y'(t) = —Bx{t) the single equation 

dy B X 
dx Ay 

This differential equation can be solved by the method of sepa- 
ration of variables: Aydy = Bxdx, so 



Ay^ = Bx^ + C, 



150- 
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Figure 3.1: Lanchester's model for the x vs. y battle. 

where C is an unknown constant. (To solve for C you must be 
given some initial conditions.) The quantity Bx^ is called the 
fighting strength of the X-men and the quantity A'lf' is called the 
fighting strength of the Y-men ("fighting strength" is not to be 
confused with "troop strength"). This relationship between the 
troop strengths is sometimes called Lanchester's square law 
and is sometimes expressed as saying the relative fight strength 
is a constant: 

Ay — Bx = constant. 



Suppose your total number of troops is some number T, where 



x{0) are initially placed in a fighting capacity and T — x{0) are 
in a support role. If your tropps outnumber the enemy then you 
want to choose the number of support units to be the smallest 
number such that the fighting effectiveness is not decreasing 
(therefore is roughly constant). The remainer should be engaged 
with the enemy in battle [M-intro]. 

A battle between three forces gives rise to the difTerential equa- 
tions 

x'{t) = -Aiy{t) - A2z{t), x{0) = xo, 
y'it) = -Bixit) - B2z{t), y(0) = 2/0, 
z'{t) = -Cix{t) - C2y{t), z(0) = zo, 

where ^4^ > 0, 5^ > 0, and > are constants and xq, and 
zq are the intial troop strengths. 
Example: Consider the battle modeled by 

x'{t) = -y{t) - z{t), x{0) = 100, 
y'{t) = -2xit)-3zit), ?/(0) = 100, 
z'{t) = -2x{t) - 3y{t), z{0) = 100. 

The Y-men and Z-men are better fighters than the X-men, in the 
sense that the coefficient of z in 2nd DE (describing their battle 
with y) is higher than that coefficient of x, and the coefficient of 
y in 3rd DE is also higher than that coefficient of x. However, 
as we will see, the worst fighter wins! 
Taking Laplace transforms, we obtain the system 

sX{s) + Y{s) + Z{s) = 100 
2X(s) + sY{s) + 3Z{s) = 100, 
2X{s) + SY{s) + sZ{s) = 100, 



which we solve by row reduction using the augmented matrix 

/ s 1 1 100 \ 



This means X{s) = ^^0^ and Y{s) = Z{s) = ^^0^. Taking 



inverse LTs, we get the solution: x{t) = 40e* + 60e~^* and y{t) = 
z{t) = — 20e* + 120e~^*. In other words, the worst fighter wins! 

In fact, the battle is over ai t = log(6)/5 = 0.35... and at this 
time, x{t) = 71.54.... Therefore, the worst fighters, the X-men, 
not only won but have lost less than 30% of their men! 

Exercise: A battle is modeled by 



(a) Write the solutions in parameteric form, (b) Who wins? 
When? State the losses for each side. 
Use SAGE to solve this. 




This has row-reduced echelon form 




x' = -4y, x{0) = 150, 
y' = -X, yiO) = 40. 




Figure 3.2: Lanchester's model for the x vs. y vs z battle. 



3.2 The Gauss elimination game and appli- 
cations to systems of DEs 



This is actually a lecture on solving systems of equations using 
the method of row reduction, but it's more fun to formulate it 
in terms of a game. 

To be specific, let's focus on a 2 x 2 system (by "2 x 2" I mean 
2 equations in the 2 unknowns x, y): 



ax -\-by = ri 
cx -\- dy = r2 



(3.1) 



Here a, 6, c, c?, ri, r2 are given constants. Putting these two equa- 
tions down together means to solve them simultaneously, not in- 
dividually. In geometric terms, you may think of each equation 
above as a line the the plane. To solve them simultaneously, 
you are to find the point of intersection (if it exists) of these two 
lines. Since a, b, c, d, ri, r2 have not been specified, it is conceiv- 
able that there are 

• no solutions (the lines are parallel but distinct), 

• infinitely many solutions (the lines are the same), 

• exactly one solution (the lines are distinct and not parallel) . 

"Usually" there is exactly one solution. Of course, you can solve 
this by simply manipulating equations since it is such a low- 
dimensional system but the object of this lecture is to show you 
a method of solution which is "scalable" to "industrial-sized" 
problems (say 1000 x 1000 or larger). 
Strategy: 

Step 1: Write down the augmented matrix oi (3.1): 



This is simply a matter of stripping off the unknowns and record- 
ing the coefficients in an array. Of course, the system must be 
written in "standard form" (all the terms with "x" get aligned 
together, ...) to do this correctly. 




Step 2: Play the Gauss elimination game (described below) to 
computing the row reduced echelon form of A, call it B say. 
Step 3: Read off the solution from the right-most column of B. 

The Gauss Elimination Game 

Legal moves: These actually apply to any mxn matrix A with 
m < n. 

1. Ri Rj-. You can swap row i with row j. 

2. cRi Ri (c 7^ 0): You can replace row i with row i mul- 
tiplied by any non-zero constant c. (Don't confuse this c 
with the c in (3.1)). 

3. cRi + Rj Ri (c 7^ 0): You can replace row i with row i 
multiplied by any non-zero constant c plus row j-, j ^ i- 

Note that move 1 simply corresponds to reordering the system 
of equations (3.1)). Likewise, move 2 simply corresponds to 
scaling equation i in (3.1)). In general, these "legal moves" 
correspond to algebraic operations you would perform on (3.1)) 
to solve it. However, there are fewer symbols to push around 
when the augmented matrix is used. 

Goal: You win the game when you can achieve the following 
situation. Your goal is to find a sequence of legal moves leading 
to a matrix B satisfying the following criteria: 

1. all rows of B have leaading non-zero term equal to 1 (the 
position where this leading term in B occurs is called a 
pivot position), 

2. B contains as many O's as possible 

3. all entries above and below a pivot position must be 0, 



4. the pivot position of the i row is to the left and above 
the pivot position of the {i + 1)** row (therefore, all entries 
below the diagonal of B are 0). 

This matrix B is unique (this is a theorem which you can find 
in any text on elementary matrix theory or linear algebra^) and 
is called the row reduced echelon form of A, sometimes written 
rref{A). 

Two comments: (1) If you are your friend both start out play- 
ing this game, it is likely your choice of legal moves will differ. 
That is to be expected. However, you must get the same result 
in the end. (2) Often if someone is to get "stuck" it is becuase 
they forget that one of the goals is to "kill as many terms as 
possible (i.e., you need B to have as many O's as possible). If 
you forget this you might create non-zero terms in the matrix 
while killing others. You should try to think of each move as 
being made in order to to kill a term. The exception is at the 
very end where you can't kill any more terms but you want to 
do row swaps to put it in diagonal form. 

Now it's time for an example. 

Example: Solve 




3 
6 



(3.2) 



The augmented matrix is 




^For example, [B-rref] or [H-rref|. 



Figure 3.3: lines x + 2y — 3, Ax + 5y — 6 in the plane. 



—4Ri + i?2 ^ i?2, which leads to 

— (l/3)i?2 R2, which leads to 

— 2i?2 + i?i ^ which leads to 

Now we are done (we won!) since this matrix satisfies all the 
goals for a eow reduced echelon form. 
The latter matrix corresponds to the system of equations 

( x + Oy = -1 . . 

\Ox + y = 2 ^"^"^^ 

Since the "legal moves" were simply matrix analogs of algebraic 
manipulations you'd appy to the system (3.2), the solution to 
(3.2) is the same as the solution to (3.3), whihc is obviously 
X = —1,2/ = 2. You can visually check this from the graph 
given above. 
To find the row reduced echelon form of 





(I 2 3\ 

V 4 5 6 y 

using SAGE , just type the following: 

SAGE 

sage: MS = MatrixSpace (QQ, 2 , 3 ) 
sage: A = MS ( [ [1, 2, 3] , [4, 5, 6] ] ) 

sage: A 
[1 2 3] 
[4 5 6] 

sage: A . echelon_f orm ( ) 

[1 0-1] 

[012] 



Solving systems using inverses 

There is another method of solving "square" systems of linear 
equations which we discuss next. 
One can rewrite the system (3.1) as a single matrix equation 

or more compactly as 

AX = f, (3.4) 

where X = ^ ^ ^ and r = ^ ^ • How do you solve (3.4)? 
The obvious this to do ( "divide by A" ) is the right idea: 



^ ) = X = A-^r. 

y J 



Here A is a matrix with the property that A A = /, the 
identity matrix (which satisfies IX = X). 

If A~^ exists (and it usually does), how do we compute it? 
There are a few ways. One, if using a formula. In the 2x2 case, 
the inverse is given by 

fa _ 1 f d -h\ 

\c d ) ad - be \ —c a J ' 

There is a similar formula for larger sized matrices but it is so 
unwieldy that is is usually not used to compute the inverse. In 
the 2x2 case, it is easy to use and we see for example. 




I f 5 -2\ f -5/3 2/3 



-3 V -4 1 y V 4/3 -1/3 
To find the inverse of 



1 2 
4 5 



using SAGE , just type the following: 

SAGE _ 



sage: MS = MatrixSpace (QQ, 2 , 2 ) 
sage: A = MS([[1,2], [4,5]]) 
sage: A 
[1 2] 
[4 5] 

sage : A' (-1) 
[-5/3 2/3] 
[ 4/3 -1/3] 



A better way to compute ^4 ^ is the following. Compute the 
row reduced echelon form of the matrix (^, /), where / is the 



identity matrix of the same size as A. This new matrix will 
be (if the inverse exists) {I,A~^). You can read off the inverse 
matrix from this. 

Here is an example. 

Example Solve 



x + 2y = 3 
Ax -\- by = 6 



using matrix inverses. 
This is 



4 5) ( y ) ( 6 



so 



{y) (45) (e)- 

To compute the inverse matrix, apply the Gauss elimination 
game to 

/ 1 2 1 \ 

V 4 5 1 y 

Using the same sequence of legal moves as in the previous ex- 
ample, we get 

/ 1 2 1 \ 

— 4i?i + i?2 — ^ i?2, which leads to 



— (l/3)i?2 — ^ -R2, which leads to y 

— 2i?2 + -Ri — ^ which leads to 
Therefore the inverse is 



yo -3 -4 1 J 
12 1 \ 
1 4/3 -1/3 J 
f 1 -5/3 2/3 
1 4/3 -1/3 



,_i ^ / -5/3 2/3 \ 

V 4/3 -1/3 ; ■ 

Now, to solve the system, compute 



fx\_/l 2y'/3\_/-5/3 2/3W3\ 

{y) [i 5) [e) I, 4/3 -1/3 ){(>) { 2 J 

To make SAGE do the above computation, just type the fohow- 
ing: 

SAGE 

sage: MS = MatrixSpace (QQ, 2 , 2 ) 
sage: A = MS ([[1,2], [4,5]]) 
sage: V = VectorSpace (QQ, 2 ) 
sage : v = V ( [ 3 , 6 ] ) 
sage : A' (-1) *v 
(-1, 2) 



Application: Solving systems of DEs 

Suppose we have a system of DEs in "standard form" 

( x' = ax + by + f{t), x{0) = xq, 
\ y' = cx + dy -\- g{t), y{0) = y^^ 

where a, 6, c, xq, yo are given constants and f{t)^g{t) are given 
"nice" functions. (Here "nice" will be left vague but basically 
we don't want these functions to annoy us with any bad be- 
haviour while trying to solve the DEs by the method of Laplace 
transforms.) 

One way to solve this system if to take Laplace transforms of 
both sides. If we let 



X{s) = C[xms),Y{s) = C[yms),Fis) = = >Cb(t)](s), 



then (3.5) becomes 

r sX{s)-xo = aX{s) + bY{s) + F{s), 

\ sY{s)-yo = cX{s) + dY{s) + Gis). ^"^'^^ 

This is now a 2 x 2 system of hnear equations in the unknowns 
X(s),Y{s) with augmented matrix 



A = 



s — a —b F{s) + Xq 
—c s — d G{s) + yo 



Example: Solve 

( x' = -y + 1, x{0) =0, 
\y' = -x + t, 1/(0) =0, 

The augmented matrix is 

/.I l/s\ 
[is III? j ■ 

The row reduced echelon form of this is 

/ 1 l/s2 

vol 

Therefore, Xi^s) = and Yi^s) = 0. Taking inverse Laplace 
transforms, we see that the solution to the system is x{i) = i 
and y{i) = 0. It is easy to check that this is indeed the solution. 

To make SAGE compute the row reduced echelon form, just 
type the following: 




SAGE _ 

sage: R = PolynomialRing (QQ, "s") 

sage: F = FractionField (R) 

sage: s = F . gen ( ) 

sage: MS = MatrixSpace (F, 2, 3) 

sage: A = MS([[s,l,l/s],[l,s,l/s"2]]) 

sage: A. echelon_f orm ( ) 

[ 1 l/s"2] 

[010] 



To make SAGE compute the Laplace transform, just type the 
following: 

SAGE 

sage: maxima ("laplace (1, t, s) ") 
1/s 

sage: maxima (" laplace (t , t, s )" ) 
l/s"2 



To make SAGE compute the inverse Laplace transform, just 
type the following: 

SAGE 

sage: maxima (" ilt ( 1/s ' 2, s, t )" ) 

t 

sage: maxima ("ilt (1/ (s^2+l) , s,t) ") 
sin (t) 



Example: Solve 



( x' = x{0) = 400, 

\y' = -X, y{Q) = 100, 



This models a battle between "ar-men" and "?/-men" , where the 
"x-men" die off at a higher rate than the "y-men" (but there 
are more of them to begin with too). 
The augmented matrix is 



A = 



s 4 400 
1 s 100 



The row reduced echelon form of this is 

1 



1 — 



100(s-4) 
s2-4 



Therefore, 



X(s) =400^^-200^ -, y(s) = 100^^-200^ -. 

Taking inverse Laplace transforms, we see that the solution to 
the system is x{t) = 400cosh(2t) - 200sinh(2t) and y{t) = 
100cosh(2t) - 200sinh(2t). The "x-men" win and, in fact, 

a:(0.275) = 346.4102..., |/(0.275) = -0.1201... . 

Question: What is x{t)^ — 4y{ty7 (Hint: It's a constant. Can 
you explain this?) 

To make SAGE plot this just type the following: 

SAGE 

sage: f = lambda x: 400*cosh (2*x) -200*sinh (2*x) 
sage: g = lambda x: 100*cosh (2*x) -200*sinh {2*x) 
sage: P = plot (f , 0, 1) 
sage: Q = plot (g, 0,1) 
sage: show(P+Q) 



sage : g ( . 2 75 ) 

-0 . 12017933629675781 
sage: f(0.275) 

346.41024490088557 




Figure 3.4: curves x{t) = 400cosh(2i) - 200sinh(2i), y{t) = 100cosh(2i) - 
200smh(2i) along the i-axis. 

Example: The displacement from equilibrium (respectively) 
for coupled springs attached to a wall on the left 

I coupled springs , 

I \/\/\/\/\— Imassll \/\/\/\/\/ |mass2| 



springl spring2 



is modeled by the system of 2nd order ODEs 

mix'l + {ki + k2)xi - k2X2 = 0, 1712X2 + k2{x2 - xi) = 0, 

where xi denotes the displacement from equilibrium of mass 
1, denoted mi, X2 denotes the displacement from equilibrium 
of mass 2, denoted m2, and /ci, k2 are the respective spring 
constants [CS-rref]. 

As another illustration of solving linear systems of equations to 
solving systems of linear 1st order DEs, we use SAGE to solve the 
above problem with mi = 2, m2 = 1, /ci = 4, k2 = 2, xi(0) = 3, 
x[{0) = 0, X2(0) = 3, x'^iO) = 0. 

Soln: Take Laplace transforms of the first DE (for simplicity 
of notation, let x = xi, y = X2): 

SAGE +MaxirrLa 

sage: _ = maxima. eval { "x2 (t) := dlff{x(t),t, 2)") 
sage : maxima ("laplace {2*x2 (t)+6*x(t)-2*y(t) ,t,s) ") 

2* (-?%at Cdiff (x(t) , t, 1) ,t=0) +s"2*?%laplace (x(t),t,s)-x(0)*s) -2*?%laplace (y (t) , t, s) +6*?%laplace (x(t) , t, s) 



This says -2x[{0) + 2s^ * Xi{s) -2sxi{0) -2X2{s) + 2Xi{s) = 
(where the Laplace transform of a lower case function is the 
upper case function). Take Laplace transforms of the second 
DE: 

SAGE+Maxima 

sage: _ = maxima. eval {"y2 (t) := diff{y(t),t, 2)") 
sage : maxima ("laplace {y2 (t ) +2*y (t) -2*x (t ) , t , s) " ) 

-?%at Cdiff (y (t) , t, 1) , t=0 ) +s "2*?%laplace (y (t) , t, s) +2*?%laplace (y (t) , t, s) -2*?%laplace (x {t ) , t, s) -y (0) *s 



This says 5^X2(5) + 2X2(5) - 2Xi(s) - 3s = 0. Solve these two 
equations: 



SAGE 



sage: s,X,Y = var('s X Y' ) 

sage: eqns = [ {2*s"2+6) *X-2*Y == 6*s, -2*X +(s"2+2)*Y == 3*s] 
sage: solve (eqns, X,Y) 

[ [X == (3*s'3 + 9*s)/(s~4 + 5*s'2 + 4), 
Y == (3*3^3 + 15*s)/(s"4 + 5*s'2 + 4)]] 



This says Xi{s) = (3s^ + 9s)/(s^ + Ss^ + 4), X2(s) = (3s^ + 
15s) /(s^ + 5s^ + 4). Take inverse Laplace transforms to get the 
answer: 

SAGE 

sage: s,t = var('s t' ) 

sage: inverse_laplace ( (3*s'~3 + 9*s)/(s"4 + 5*s"2 + 4),s,t) 
cos{2*t) + 2*cos(t) 

sage: inverse_laplace ( (3*s"3 + 15*s)/(s'4 + 5*s"2 + 4),s,t) 
4*cos(t) - cos(2*t) 



Therefore, xi{t) = cos(2t) + 2cos(t), X2{t) = 4cos(t) - cos(2t). 
Using SAGE , this can be plotted parametrically using 



. SAGE 



sage: P = parametric_plot { [cos (2*t) + 2*cos (t) , 4*cos (t) - cos (2*t) ] , 0, 3) 
sage: show(P) 



You can also try 



- SAGE +Maxima - 



sage . : maxima .plot 2d ('cos{2*x) + 2*cos{x)','[x,0,l]',' [plot_f ormat , openmath] ' ) 



for the output of a slightly different looking plotting program. 



-1 



1 



2 



3 




-1.- 



Figure 3.5: curves x{t) = cos{2 *t) + 2* cos{t), y{t) = 4 * cos{t) — cos{2 * t) 
along the t-axis. 

Exercise: Solve 

X + 2y + z = 1 
-x + 2y-z = 2 
y^2z = 3 

using (a) row reduction and SAGE, (b) matrix inverses and 
SAGE. 



3.3 Eigenvalue method for systems of DEs 



Motivation 

First, we sliall try to motivate the study of eigenvalues and 
eigenvectors. This section hopefully will convince you that 

• diagonal matrices are wonderful, 

• conjugation is very natural, 

• if our goal in life is to conjugate a given square matrix ma- 
trix into a diagonal one, then eigenvalues and eigenvectors 
are also natural. 

Diagonal matrices are wonderful: We'll focus for simplicity on 
the 2x2 case, but everything applies to the general case. 

• Addition is easy: 

/^aiO\ /6iO^/ai + 6i \ 
\ a2 J \0 b2 J V 02 + 62; 

• Multiplication is easy: 

ai \ ^ bi \ _ / ai-bi \ 
a2 J ' \ b2 J ~ \ a2 • 62 7 ■ 

• Powers are easy: 

\0 a2) ~V0 ^2)' 



• You can even exponentiate them: 




So, diagonal matrices are wonderful. 



Conjugation is natural. You and your friend are piloting a rocket 
in space. You handle the controls, your friend handles the map. 
To communicate, you have to "change coordinates" . Your coor- 
dinates are those of the rocketship (straight ahead is one direc- 
tion, to the right is another). Your friends coordinates are those 
of the map (north and east are map directions). Changing co- 
ordinates corresponds algebraically to conjugating by a suitable 
matrix. Using an example, we'll see how this arises in a specific 
case. 

Your basis vectors are 

i;i = (l,0), V2 = (0,l), 

which we call the "^'-space coordinates", and the map's basis 
vectors are 



w;i = (l,l), w;2 = (l,-l), 
which we call the "lu-space coordinates" . 




Figure 3.6: basis vectors vi,V2 and wi,W2. 



For example, the point (7, 3) is, in -y-space coordinates of course 
(7, 3) but in the w-space coordinates, (5, 2) since + 2w2 = 

7vi + 3v2- Indeed, the matrix A = ( ] \ | sends I ^ I to 



1 -1 



7 

3 



Suppose we flip about the 45^ line (the "diagonal") in each 
coordinate system. In the 'U-space: 

avi + bv2 I — > bvi + av2, 

(:)-(;;)(:)■ 



In other words, in 'U-space, the "flip map" is 
In the w;-space: 

wvi + WV2 I — > awi — hw2-, 

(;)-(; -°.)(:) 



1 

1 



In other words, in ly-space, the "flip map" is 



1 « V 

Conjugating by the matrix A converts the "flip map" in w- 
space to the the "flip map" in 'U-space: 



Eigenvalues are natural too 

Given a matrix A, is there a basis of the underlying space in 
which the matrix is diagonal? Given how "wonderful" diagonal 
matrices are, it seems clear we should flnd this basis and these 
diagonal entries. 

Fact When the diagonal entries are distinct, the basis elements 
are the eigenvectors and the diagonal elements are the eigenval- 
ues. 



Since this section is only intended to be motivation, we shall 
not prove this here (see any text on linear algebra, for example 
[B-rref] or [H-rref]). 

SAGE 

sage: MS = MatrixSpace (CC, 2 , 2 ) 
sage: A = MS ( [ [0, 1] , [1, 0] ] ) 
sage: A. eigenspaces ( ) 

[ 

(1.00000000000000, [ 

(1.00000000000000, 1 .00000000000000) 

] ) , 

(-1.00000000000000, [ 

(1. 00000000000000, -1 . 00000000000000) 

] ) 

] 



Solution strategy 



PROBLEM: Solve 



x' = ax -\- by, x{0) = xq, 
y' = cx + dy, y{0) = yo. 



soln: Let 



A = 



a b \ 
^ c d J 

In matrix notation, the system of DEs becomes 



X' = AX, X{0)=(^2)' 



where X = X(t) = ( ^[^} ] ■ In a similar manner to how we 

V J 

solved homogeneous constant coefficient 2nd order ODEs ax" + 
bx' + ex = by using "Euler's guess" x = Ce*^*, we try to guess 
an exponential: X(t) = ce^* (A is used instead of r to stick with 
notational convention; c in place of C since we need a constant 
vector). Plugging this guess into the matrix DE X' = AX gives 
Ace^* = Ace^^, or (cancelling e^*) 

Ac = Ac. 

This means that A is an eigenvalue of A with eigenvector c. 

• Find the eigenvalues. These are the roots of the character- 
istic polynomial 

p{X) = det ^ ^ ~ = - {a + d)X + {ad - be). 

Call them Ai, A2 (in any order you like). 



You can use the quadratic formula, for example to get them: 



a + d {a -\- dy — 4(a(i — be) a + d y^(a + d)"^ — 4:{ad — be) 

^1 = — 7^ 1 7^ 5 ^^2 = — 7^ 



• Find the eigenvectors. If 6 7^ then you can use the for- 
mulas 



\ Xi-a J ' V A2 - a y ' 

In general, you can get them by solving the eigenvector 
equation Av = Xv. 

SAGE 

sage: MS = MatrixSpace (CC, 2, 2) 
sage: A = MS{ [ [1,2] , [3,4] ] ) 
sage: A. eigenspaces ( ) 



(-0.372281323269014, [ 

(1. 00000000000000, -0. 45742710775 6338) 
] ) , 

(5.37228132326901, [ 

(1.00000000000000, 1. 457 42710775 634) 

] ) 

] 



• Plug these into the following formulas: 

(a) Ai 7^ A2, real: 

( yft) ) ^ ciiTi exp(Ait) + e2V2Qw{^2t). 

(b) Ai = A2 = A, real: 

( yft) ) ^ ^^^^ exp(At) + e2{vit + p) exp(At), 



where p is any non-zero vector satisfying [A — XI)p = 

Vi. 

(c) Ai = Q! + i/3, complex: write vi = ui -\- iu2, where ui 
and U2 are both real vectors. 




= ci[ex.p{at) cos{(3t)ui — ex.p{at) sm{f3t)u2\ 



+C2[— exp{at) cos{/3t)u2 — exp{at) sin{/3t)ui\. 



Examples 

Example 3.3.1. Solve 

x'{t))=x{t)-y{t), y'{t)=4x{t)^y{t), x{0) = -1, y{0) = l. 
Let 

and so the characteristc polynomial is 

p{x) = det(74 — xl) = — 2x + 5. 
The eigenvalues are 

Ai = l + 2i, A2 = l-2i, 

so a = I and (3 = 2. Eigenvectors vi,V2 are given by 

though we actually only need to know vi . The real and imaginary 
parts of vi are 



The solution is then 

/ x{t) \ _ f — ci exp(t) cos(2t) + C2 exp(t) sin(2t) \ 
\^ ylt) y V -2ciexp(t)sin(2t) - 2c2 exp(t) cos(2t), J 

so x{t) = —c\ exp(t) cos(2t)+C2exp(t) sin(2t) andy{t) = — 2ci exp(t) sin(2t) — 
2c2exp(t) cos(2t). 

Since x{0) = —1, we solve to get ci = 1. Since y{0) = 1, 
we get C2 = —1/2. The solution is: x{t) = — exp(t) cos(2t) — 
I exp(t) sin(2t) and y{t) = — 2exp(t) sin(2t) + exp(t) cos(2t). 

Example 3.3.2. Solve 

x\t) = -2x{t) + Syit), y\t) = -?>x{t) + Ay{t). 

Let 

and so the characteristc polynomial is 

p{x) = det(A - xl) = x^ - 2x + 1. 

The eigenvalues are 

Ai = A2 = 1. 
An eigenvector vi is given by 

•-.=(?) 

Since we can multiply any eigenvector by a non-zero scalar and 
get another eigenvector, we shall use instead 



Letp ~ \^gj non-zero vector satisfying (A — XI)p = vi. 

This means 

There are infinitely many possibly solutions but we simply take 
r = Oands = 1/3, so 



p = 





1/3 



The solution is 



(^(tl)=^'(i)<=^pw+'=4i)*+(i;3))-p«' 

orxit) = ci exp(t)+C2t exp(t) andyit) = ci ex.p{t) -\- ex.p{t) -\- 
C2t exp(t) . 



Exercises: Use SAGE to find eigenvalues and eigenvectors of 
both 

and 



3.4 Electrical networks using Laplace trans- 
forms 

Suppose we have an electrical network (i.e., a series of electri- 
cal circuits) involving emfs (electromotive forces or batteries), 
resistors, capacitors and inductors. We use the following "dic- 
tionary" to translate between the diagram and the DEs. 



EE object 


term in DE 
(the voltage drop) 


units 


symbol 


charge 


q = J i{t) dt 


coulombs 




current 


i = q' 


amps 




emf 


e = e{t) 


volts V 




resistor 


Rq' = Ri 


ohms Q 




capacitor 


C~^q 


farads 


— 1 1— 


inductor 


Lq" = Li' 


henries 





Kirchoff's First Law. The algebraic sum of the currents trav- 
elling into any node is zero. 

KirchojJ's Second Law: The algebraic sum of the voltage drops 
around any closed loop is zero. 

Example 1: Consider the simple RC circuit given by the fol- 
lowing diagram. 

According to Kirchoff's 2"^^ Law and the above "dictionary", 
this circuit corresponds to the DE 

g' + 5g = 2. 

The general solution to this is q{t) = 1 + ce~^*, where c is a 
constant which depends on the initial charge on the capacitor. 



1 ohm 



T 



2 volts 



1/5 farads 
Figure 3.7: A simple circuit. 



□ 



Aside: The convention of assuming that electricity flows from 
positive to negative on the terminals of a battery is referred 
to as "conventional flow". The physically-correct but opposite 
assumption is referred to as "electron flow". We shall assume 
the "electron flow" convention. 

Example 2: Consider the network given by the following di- 
agram. 

1 ohm 

> — ^AA >i 



T 

nodel ^ 



2 volts 



1/5 farads 
—I I ^ 



node 2 



5 ohms 



Figure 3.8: A network. 



Assume the initial charges are 0. 



One difference between this circuit and the one above is that 
the charges on the three paths between the two nodes (labeled 
node 1 and node 2 for convenience) must be labeled. The charge 
passing through the 5 ohm resistor we label qi, the charge on 
the capacitor we denote by q2, and the charge passing through 
the 1 ohm resistor we label q^. 

There are three closed loops in the above diagram: the "top 
loop", the "bottom loop", and the "big loop". The loops will 
be traversed in the "clockwise" direction. Note the "top loop" 
looks like the simple circuit given in Example 1 but it cannot 
be solved in the same way, since the current passing through 
the 5 ohm resistor will affect the charge on the capacitor. This 
current is not present in the circuit of Example 1 but it does 
occur in the network above. 

Kirchoff's Laws and the above "dictionary" give 

9^ + 5^2 = 2, gi(0) = 0, 
5^1-5^2 = 0, 92(0) = 0, 
5gi + gJ = 2, ^3(0) = 0. 

Notice the minus sign in front of the term associated to the 
capacitor (— 5^2)- This is because we are going clockwise, against 
the "direction of the current" . Kirchoff's 1** law says q'^ = q'i+q2- 
Since qi{0) = ^2(0) = ^3(0) = 0, this implies ^3 = ^1 + q2- After 
taking Laplace transforms of the 3 differential equations above, 
we get 

sQsis) + 5^2(5) = 2/s, 5sQi{s) - 5^2(5) = 0. 

Note you don't need to take th eLT of the S^'^ equation since 
it is the sum of the first two equations. The LT of the above 
+ 52 = qs (Kirchoff's law) gives Qi{s) + ^2(5) - Qsis) = 0. 



We therefore have this matrix equation 



Qi{s) 

Q2{S] 

•1/ \Q^{s) 

The augmented matrix describing this system is 






The row-reduced echelon form is 



/ 1 2/(s3 + 6s2) 
1 2/(s2 + 6s) 
\0 1 2(s + l)/(s3 + 6s2) 



Therefore 



Qi{s) = 



This imphes 



Q2[S) = 



2(5+1) 

s2(s + 6)' 



qi{t) = -l/18+e-^7l8+i/3, q2{t) = l/3-e-^V3, qs{t) = q2{t)+qi{t). 
□ 

This computation can be done in SAGE as well: 

SAGE 

sage : s = var ( " s " ) 

sage: MS = MatrixSpace (SymbolicExpressionRing ( ) , 3, 4) 
sage: A = MS { [ [0, 5, s, 2/s] , [5*s, 0, s, 2/s] , [1, 1, -1, 0] ] ) 
sage: B = A. echelon_f orm ( ) ; B 



[ 



2/(5*s'2) - (-2/{5*s) - 2/{5*s'~2) 
2/(5*s) - (-2/(5*3) - 2/(5*s~2))*s 
(-2/ (5*s) - 2/ (5*s"2) ) 



sage: Ql = B [0, 3] 

sage : t = var ( "t " ) 

sage: Ql . inverse_laplace (s, t) 

e^ (- (6*t) ) /18 + t/3 - 1/18 

sage : Q2 = B [1, 3] 

sage: Q2 . inverse_laplace (s, t) 

1/3 - e" (-(6*t) ) /3 

sage : Q3 = B [2, 3] 

sage: Q3 . inverse_laplace (s, t) 

-5*e" (- (6*t) ) /IB + t/3 + 5/18 



/(5*(-s/5 - 6/5) 
(5*(-s/5 - 6/5)) 
(-8/5 - 6/5) 



Example 3: Consider the network given by the following dia- 
gram. 




-> 1 



0.1 henry 



Figure 3.9: Another network. 



Assume the initial charges are 0. 



Using Kirchoff 's Laws, you get a system 

k-'i2-h = 0, 
2zi + Z2 + (0.2)z; = 6, 

(0.1)Z'3 - Z2 = 0. 

Take LTs of these three DEs. You get a 3 x 3 system in the 
unknowns I\{s) = £[ii(t)](s), /2(s) = £[i2(t)](s), and /3(s) = 
£[i3(t)](s). The augmented matrix of this system is 

1 -1-10 
2 + s/5 1 6/s 
-1 s/10 

(Check this yourself!) The row-reduced echelon form is 

I ^ ^ ^ s(s2+25s+100) * 
10 



s2+25s+100 

n n 1 300 I 

\ ^ ^ s(s2+25s+100) / 



Therefore 



1 2 3 
iJs) = h-. 

This implies 



s + 20 s + 5' 



his) = 



i^it) = 3-2e-^*-e-'"*, i2{t) = 2e-^*-2e-'"^ Z3(^) = 3-4e-^*+e 
□ 

Exercise: Use SAGE to solve for ^2(^)5 and i^{t) in the 

above problem. 



Chapter 4 



Introduction to partial 
differential equations 

4.1 Introduction to separation of variables 

A partial differential equation (PDE) is an equation satisfied by 
an unknown function (called the dependent variable) and its 
partial derivatives. The variables you differentiate with respect 
to are called the independent variables. If there is only one 
independent variable then it is called an ordinary differential 
equation. 
Examples include 

• the Laplace equation ^+|^ = 0, where u is the dependent 
variable and x,y are the independent variables, 

• the heat equation Ut = auxxi 

• and the wave equation Uu = c^Uxx- 

All these PDEs are of second order (you have to differentiate 
twice to express the equation). Here, we consider a first order 
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PDE which arises in applications and use it to introduce the 
method of solution called separation of variables. 



The transport or advection equation 

Advection is the transport of a some conserved scalar quantity 
in a vector field. A good example is the transport of pollutants 
or silt in a river (the motion of the water carries these impurities 
downstream) or traffic flow. 

The advection equation is the PDE governing the motion of 
a conserved quantity as it is advected by a given velocity field. 
The advection equation expressed mathematically is: 

_ + V • (wa) = 

where V- is the divergence and a is the velocity field of the fluid. 
Frequently, it is assumed that V • a = (this is expressed by 
saying that the velocity field is solenoidal). In this case, the 
above equation reduces to 

du „ ^ 
— + a • Vw = 0. 

ot 

Assume we have horizontal pipe in which water is flowing at a 
constant rate c in the positive x direction. Add some salt to this 
water and let u(x, t) denote the concentration (say in lbs/gallon) 
at time t. Note that the amount of salt in an interval / of the 
pipe is Jju{x,t) dx. This concentration satisfies the transport 
(or advection) equation: 

Ut + CUx = 0. 

(For a derivation of this, see for example Strauss [S-pde], §1.3.) 
How do we solve this? 



Solution 1: D'Alembert noticed that the directional derivative 
of u(x, t) in the direction v = ^^^^ (c, 1) is D^{u) = '^/^^{cUx + 
Ut) = 0. Therefore, u{x, t) is constant along the lines in the 
direction of v, and so u{x,t) = f{x — ct), for some function /. 
We will not use this method of solution in the example below 
but it does help visualize the shape of the solution. For instance, 
imagine the plot of z = f{x — ct) in {x, t, z) space. The contour 
lying above the line x = ct -\- k {k fixed) is the line of constant 
height z = f{k). □ 

Solution 2: The method of separation of variables indicates 
that we start by assuming that u{x, t) can be factored: 

u{x,t) = X{x)T{t), 

for some (unknown) functions X and T. (One can shall work 
on removing this assumption later. This assumption "works" 
because partial diflterentiation of functions like xH^ is so much 
simpler that partial differentiation of "mixed" functions like 
sin(a:^ + t^).) Substituting this into the PDE gives 

X{x)r{t) + cX\x)T{t) = 0. 

Now separate all the x's on one side and the t's on the other 
(divide by X(a;)T(i)): 

T'jt) _ _ X'jx) 
~ ~^X{x) ' 

(This same "trick" works when you apply the separation of vari- 
ables method to other linear PDE's, such as the heat equation 
or wave equation, as we will see in later lessons.) It is impossi- 
ble for a function of an independent variable x to be identically 
equal to a function of an independent variable t unless both are 



T'(t) 

constant. (Indeed, try taking the partial derivative of with 
respect to x. You get since it doesn't depend on x. Therefore, 
the partial derivative of —Cj^ is akso 0, so is a constant!) 

Therefore, = — = K, for some (unknown) constant 
K. So, we have two ODEs: 

T'(t) X'ix) 

Therefore, we can converted the PDE into two ODEs. Solving, 
we get 

T{t) = cie^\ X{x) = C2e-^^/^ 
so, u{x,t) = j[^Kt-Kx/c _ ^e-f(a;-ct)^ fQj. some constants K 
and A (where A is shorthand for C1C2; in terms of D'Alembert's 
solution, f{y) = Ae^^^y^^). The "general solution" is a sum of 
these (for various A's and K's). □ 

This can also be done in SAGE : 

SAGE 



sage: t = var("t"; 
sage: T = function { "T" , t) 
sage: K = var("K") 
sage: TO = var{"TO") 
sage: maxima . de_solve (' diff (T, t) =\ 

K*T', ['t','T'], [0,T0]) 

T=%e'^ (t*K) *T0 

sage: x = var("x") 

sage: X = f unction ( "X" , x) 

sage: c = var{"c") 

sage: XO = var("XO") 

sage: maxima . de_solve (' diff (X, x) =\ 

-c' (-1) *K*X' , ['x','X'], [0,X0] 
X=%e"- (x*K/c) *X0 

sage: solnX = maxima . de_solve (' diff (X, x) =\ 



-c" (-1) *K*X' , ['x','X'], [0,X0]) 

sage: solnX.rhsO 
%e"- (x*K/c) *X0 

sage: solnT = maxima . de_solve (' diff (T, t ) =\ 

K*T', ['t','T'], [0,T0]) 
sage: solnT.rhsO 
%e'^ (t*K) *T0 

sage: solnT . rhs ( ) *solnX . rhs ( ) 
%e' (t*K-x*K/c) *TO*XO 



Example: Assume water is flowing along a horizontal pipe at 
3 gal/min in the x direction and that there is an initial con- 
centration of salt distributed in the water with concentration of 
u{x^ 0) = e~^. Using separation of variables, find the concentra- 
tion at time t. Plot this for various values of t. 

Solution: The method of separation of variables gives the "sep- 
arated form" of the solution to the transport PDE as u{x^ t) = 
j^^Kt~Kx/c^ where c = 3. The initial condition implies 

e-^ = u{x, 0) = Ae^-o-^^/^ = Ae-^^/^ 

so ^4 = 1 and K = 3. Therefore, u{x^t) = e^^^-'^'. In other words, 
the salt concentration is increasing in time. This makes sense if 
you think about it this way: "freeze" the water motion at time 
t = 0. There is a lot of salt at the beginning of the pipe and 
less and less salt as you move along the pipe. Now go down the 
pipe in the x-direction some amount where you can barely tell 
there is any salt in the water. Now "unfreeze" the water motion. 
Looking along the pipe, you see the concentration is increasing 
since the saltier water is now moving toward you. 
This is produced using either the Maxima command 



Figure 4.1: Transport with velocity c = 3. 



Maxima 

(%il) plot3d (exp (3*t-x) , [x,0,2], [t,0,2], [grid, 12, 12] ) ; 



or the SAG E command 

SAGE 

sage: maxima .plot3d ( "exp (3*t-x) " , "[x,0,2]", "[t,0,2]",\ 
" [grid, 12 , 12 ] " , ' [plot_format, openmath] ' ) 



In both cases, wish and tcl/tk must also be installed. 
□ 

What if the initial concentration was not u{x,0) = but 
instead u{x, 0) = e^^ + 3e~^^? How does the solution to 

w^ + 3w:, = 0, w(a:,0) = e~^ + 3e-^^ (4.1) 

differ from the method of solution used above? In this case, we 
must use the fact that (by superposition) "the general solution" 
is of the form 



u{x, t) = Aie^'^'-^'l'^^ + ^26^'^^*-^/^^ + ylge^^^^-^/^) + ... , (4.2) 

for some constants Ai, Ki, .... To solve this PDE (4.1), we must 
answer the following questions: (1) How many terms from (4.2) 
are needed? (2) What are the constants Ai,Ki,...7 There are 
two terms in u{x, 0), so we can hope that we only need to use 
two terms and solve 

e-^ + 3e-^^ = u{x, 0) = Aie^^^o-^/^^ + yl2e^2(o-x/3) 

for Ai, Ki, A2, K2. Indeed, this is possible to solve: Ai = 1, 
Ki = 3, ^2 = 3, Ki = 15. This gives 

w(x,t) = e^(*-^/3) + 3ei5(*-^/^). 

Exercise: Using SAGE, solve and plot the solution to the fol- 
lowing problem. Assume water is flowing along a horizontal pipe 
at 3 gal/min in the x direction and that there is an initial con- 
centration of salt distributed in the water with concentration of 

u{x, 0) = e^. 



4.2 Fourier series, sine series, cosine series 



History: Fourier series were discovered by J. Fourier, a French- 
man who was a mathematician among other things. In fact, 
Fourier was Napolean's scientific advisor during France's inva- 
sion of Egypt in the late 1800's. When Napolean returned to 
France, he "elected" (i.e., appointed) Fourier to be a Prefect 
- basically an important administrative post where he oversaw 
some large construction projects, such as highway constructions. 
It was during this time when Fourier worked on the theory of 
heat on the side. His solution to the heat equation is basically 
what undergraduates often learn in a DEs with EVPs class. The 
exception being that our understanding of Fourier series now is 
much better than what was known in the early 1800's and some 
of these facts, like Dirichlet's theorem, are covered as well. 

Motivation: Fourier series, since series, and cosine series are 
all expansions for a function /(x), much in the same way that a 
Taylor series ao + ai{x — xq) -\- a2{x — xq)'^ + ... is an expansion. 
Both Fourier and Taylor series can be used to approximate f{x). 
There are at least three important differences between the two 
types of series. (1) For a function to have a Taylor series it 
must be differentiable^ , whereas for a Fourier series it does not 
even have to be continuous. (2) Another difference is that the 
Taylor series is typically not periodic (though it can be in some 
cases), whereas a Fourier series is always periodic. (3) Finally, 
the Taylor series (when it converges) always converges to the 
function f{x), but the Fourier series may not (see Dirichlet's 



^Remember the formula for the n-th Taylor series coefficient centered at a; = a;o 



theorem below for a more precise description of what happens). 

Definitions: Let f{x) be a function defined on an interval of 
the real line. We allow f{x) to be discontinuous but the points 
in this interval where f{x) is discontinuous must be finite in 
number and must be jump discontinuities. 

First, we discuss Fourier series. To have a Fourier series you 
must be given two things: (1) a "period" P = 2L, (2) a function 
f{x) defined on an interval of length 2L, usually we take —L < 
X < L (but sometimes < x < 2L is used instead). The Fourier 
series of f{x) with period 2L is 



00 



r/ N O.0 v^r .riTTx. , . .mix.. 
/(^) ~ y + cos(— ) + hn sm(— )J 

n=l 

where an and hn are given by the formulas^, 



1 TL7TX 

fln = y f{x) cos(-^) dx, (4.3) 



and 



1 TLTTX 

&n = y ^ fix) sin(-^) dx. (4.4) 

Next, we discuss cosine series. To have a cosine series you must 
be given two things: (1) a "period" P = 2L, (2) a function f{x) 
defined on the interval of length L, < x < L. The cosine 
series of f{x) with period 2L is 

00 

f{x) ~ y + 2^anCos(— ), 

n=l 

where a„ is given by 



^Thcsc formulas were not known to Fourier. To compute the Fourier coefficients a„, 6„ 
he used sometimes ingenious round-about methods using large systems of equations. 



2 



/ cos(— — dx. 

Jo Li 

(This formula is not in your USNA Math Tables.) The cosine 
series of f{x) is exactly the same as the Fourier series of the 
even extension of /(x), defined by 



/(x), < x < L, 
f{-x), -L <x <0. 

Finally, we define sine series. To have a sine series you must 
be given two things: (1) a "period" P = 2L, (2) a function f{x) 
defined on the interval of length L, < x < L. The sine series 
of f{x) with period 2L is 



00 

,n7TX. 



n—\ 



where 6„ is given by 

2 



K = 



L 



sin(^)/(a;) dx. 



The sine series of f{x) is exactly the same as the Fourier series 
of the odd extension of f{x), defined by 



foddix) = 



f{x), <x < L, 
-/(-x), -L <x <0. 

One last definition: the symbol ~ is used above instead of = 
because of the fact that was pointed out above: the Fourier 
series may not converge to f{x). Do you remember right-hand 
and left-hand limits from calculus 1? Recall they are denoted 
= lime^o,e>o /(a; + e) and /(x-) = lime^o.oo /(a^ - e), 



resp.. The meaning of ~ is that the series does necessarily not 
converge to the value of f{x) at every point^. The convergence 
proprties are given by the theorem below. 

Dirichlet's theorem^: Let f{x) be a function as above and 
let —L < X < L. The Fourier series of f{x), 

00 

fix) ~ y + 2^[anCos(-^) + 6nsm(-^)j, 

n=l 

(where and 6„ are as in the formulas (4.3), (4.4)) converges 
to 

fix+) + fix-) 
2 

In other words, the Fourier series of f{x) converges to f{x) only 
if f{x) is continuous at x. If f{x) is not continuous at x then 
then Fourier series of f{x) converges to the "midpoint of the 
jump" . 

Examples: (1) If f{x) = 2-\-x, —2<x<2, then the definition 
of L implies L = 2. Without even computing the Fourier series, 
we can evaluate it using Dirichlet's theorem. 

Question: Using periodicity and Dirichlet's theorem, find the 
value that the Fourier series of f{x) converges to at a: = 1, 2, 3. 

(Ans: f{x) is continuous at 1, so the FS at a; = 1 converges to /(I) = 3 by 
Dirichlet's theorem. f{x) is not defined at 2. It's FS is periodic with period 4, so at 
X = 2 the FS converges to /(2+)+/(2-) ^ o±4 ^ 2. f{x) is not defined at 3. It's FS 
is periodic with period 4, so at a; = 3 the FS converges to = 1±1 = 1.) 

The formulas (4.3) and (4.4) enable us to compute the Fourier 
series coefficients ao, fln and 6„. (We skip the details.) These 

^Fourier believed his series converged to the function in the early 1800's but we now 

know this is not always true. 
^Pronounced "Dear-ish-lay" . 



formulas give that the Fourier series of f{x) is 



A °° 




n=l 



The Fourier series approximations to /(x) are 



4 ,7TX, sinf^TTx) 
So = 2, Si = 2+- sin(— ), ^2 = 2+4 ^-2 

TT 2 TT 



sin (tt x) 



TT 



The graphs of each of these functions get closer and closer to 
the graph of f{x) on the interval — 2 < x < 2. For instance, the 
graph of f{x) and of 5*8 are given below: 

Notice that f{x) is only defined from — 2 < a; < 2 yet the Fourier 
series is not only defined everywhere but is periodic with period 
P = 2L = 4. Also, notice that 5*8 is not a bad approximation to 



This can also be done in SAGE . First, we define the function. 

SAGE 

sage: f = lambda x:x+2 

sage: f = Piecewise ( [ [ (-2, 2) , f ] ] ) 



This can be plotted using the command f . plot ( ) . show ( ) . Next, 
we compute the Fourier series coefficients: 

SAGE 

sage: f . fourier_series_cosine_coef f icient (0, 2) # a_0 
4 

sage: f. fourier_series_cosine_coeff icient (1, 2) # a_l 



sage: f . f ourier_series_cosine_coef f icient (2 , 2 ) # a_2 



sage: f . fourier_series_cosine_coeff icient (3, ) # a_3 



X 

Figure 4.2: Graph of f{x) and a Fourier series approximation of f{x). 





sage: f . f ourier_series_sine_coef f icient (1, 2) # b_l 
4/pi 

sage: f . fourier_series_sine_coef ficient (2, ) # b_2 

-2/pi 

sage: f . fourier_series_sine_coef ficient (3, 2) # b_3 
4/ (3*pi) 



Finally, the partial Fourier series and it's plot verses the function 
can be computed using the following SAGE commands. 

I 1 

sage: f . fourier_series_partial_sum (3, 2) 
-2*sin (pi*x) /pi + 4*sin (pl*x/2) /pi + 2 

sage: PI = f .plot_fourier_series_partial_sum (15, 2, -5, 5, linestyle=" : " ) 



The plot (which takes 15 terms of the Fourier series) is given 
below. 




■1- 



Figure 4.3: Graph of f{x) = x + 2 and a Fourier series approximation, L = 2. 



(2) This time, let's consider an example of a cosine series. In 
this case, we take the piecewise constant function f{x) defined 
on < x < 3 by 



1, 0<x<2, 
-1, 2 < X < 3. 



We see therefore L = 3. The formula above for the cosine series 
coefficients gives that 



n=l 




UTTX 



The first few partial sums are 



^2 = 1/3 + 2 



a/3 cos (I TTX) 



TT 



^3 = 1/3 + 2 



\/3cos(|7rx) -\/3 cos (I TT x) 



TT 



TT 



As before, the more terms in the cosine series we take, the better 
the approximation is, for < x < 3. Comparing the picture 
below with the picture above, note that even with more terms, 
this approximation is not as good as the previous example. The 
precise reason for this is rather technical but basically boils down 
to the following: roughly speaking, the more differentiable the 
function is, the faster the Fourier series converges (and therefore 
the better the partial sums of the Fourier series will approximate 
f{x)). Also, notice that the cosine series approximation is an 
even function but f{x) is not (it's only defined from < x < 3). 
For instance, the graph of f{x) and of are given below: 

(3) Finally, let's consider an example of a sine series. In this 
case, we take the piecewise constant function f{x) defined on 
< x < 3 by the same expression we used in the cosine series 
example above. 

Question: Using periodicity and Dirichlet's theorem, find the 
value that the sine series of f{x) converges to at a: = 1,2,3. 

(Ans: f{x) is continuous at 1, so the FS at x = 1 converges to /(I) = 1. f{x) is 



Figure 4.4: Graph of f{x) and a cosine series approximation of f{x). 



not continuous at 2, so at x = 2 the SS converges to ) = /( 2+)+f{2 ) 

-1+1 



0. f{x) is not defined at 3. It's SS is periodic with period 6, so at x = 3 
the SS converges to /°dd(3-)+/odd(3+) ^ -^fl ^ 

The formula above for the sine series coefficients give that 

, cos (nTr) — 2 cos (I nvr) + 1 mttx, 
fi^) = E 2 -^-^ sm — - . 

n=l 

The partial sums are 

sin(l/3 7rx) sin(|7rx) 



^2 = 2 '-^ + 3 

TT TT 

^sinf^TTx) ^sin(|7rx) ^, sin(7ra;) 
^3 = 2 ^ ^ + 3 ^ ^ - 4/3 ^ ... 

TT TT TT 

These partial sums iS^, as n ^ oo, converge to their limit about 
as fast as those in the previous example. Instead of taking only 
10 terms, this time we take 40. Observe from the graph below 
that the value of the sine series at x = 2 does seem to be ap- 
proaching 0, as Dirichlet's Theorem predicts. The graph of f{x) 
with is 
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Figure 4.5: Graph of f{x) and a sine series approximation of f{x). 



Exercise: Let f{x) = x^, —2 < x <2 and L = 2. Use SAGE to 
compute the first 10 terms of the Fourier series, and plot the 
corresponding partial sum. Next plot the partial sum of the 
first 50 terms and compare them. 

Exercise: What mathematical results do the following SAGE commands 
give you? In other words, if you can seen someone typing these 
commands into a computer, explain what problem they were 
trying to solve. 

I SAGE 

sage : x = var ( "x" ) 

sage : f = lambda x : 0*x 

sage : f 1 = lambda x : -x" 

sage : f 2 = lambda x : x"0 

sage: f = Piecewise ( [ [ (-2, 0) , f 1 ] , [ (0, 3/2) , f 0] , [ {3/2, 2) , f 2 j ] ) 

sage: PI = f .plot (} 

sage: alO = [f . f ourier_series_cosine_coef f icient (n, 2) for n in range (10) ] 

sage: blO = [f . fourier_series_sine_coeff icient {n, 2) for n in range (10) ] 

sage: fslO = alO[0]/2 + sum( [alO [i] *cos (i*pi*x/2) for i in 
range(l,10)]) + sum{ [blO [i] *sin (i+pi*x/2) for i in range(lO)]) 

sage: P2 = f slO .plot (-4, 4, linestyle=" : " ) 

sage: (P1+P2) .show() 
sage : 

sage: a50 = [f . fourier_series_cosine_coeff icient {n, 2 ) for n in range(50)] 

sage: b50 = [f . fourier_series_sine_coefficient (n, 2) for n in range (50)] 

sage: fs50 = a50[0]/2 + sum( [a50 [i] *cos {i*pi*x/2) for i in 
range (1,50)]) + sum ( [b50 [ i ] *sin (i*pi*x/2 ) for i in range (50 ) ] ) 

sage: P3 = f s50. plot (-4 , 4 , linestyle=" — ") 

sage: (P1+P2+P3) . show () 

sage : alOO = [ f . fourier_series_cosine_coeff icient (n, 2 ) for n in range { 100 ) ] 

sage : blOO = [ f . fourier_series_sine_coeff icient (n, 2 ) for n in range (100 ) ] 

sage: fslOO = al00[0]/2 + sum( [alOO [i] *cos (i*pi*x/2) for i in 

range (1, 100) ] ) + sum ( [bl 00 [ i ] *sin ( i *pi*x/2 ) for i in range(lOO)]) 

sage: P3 = f slOO .plot (-4, 4, linestYle=" — ") 

sage: (P1+P2+P3) .show{) 

sage : 



4.3 The heat equation 



The deep study of nature is the most fruitful source 
of mathematical discoveries. 

- Jean- Baptist- Joseph Fourier 



The heat equation with zero ends boundary conditions models 
the temperature of an (insulated) wire of length L: 



Here u{x, t) denotes the temperature at a point x on the wire at 
time t. The initial temperature f{x) is specified by the equation 




u{x,0) = f{x). 



Method: 



• Find the sine series of f{x): 



oo 



/W~EM/)sm(^) 



n=l 



• The solution is 



oo 



U 



Example: Let 



/(^) = 



Then L = n and 



-1, 0<x<7r/2, 
2, 7r/2 < X <7r. 



2 r ^ ^ , , 2 cos(/77r) - 3 cos(f?77r) + 1 

^n(/) = — fix) sm(nx)dx = —2 

TT Jo 

Thus 



2 6 2 

/(x) ~ bi{f) sm{x)+b2{f) sin(2a;)+... = — sin(a;) — sin(2a;)+— sin(3a; 

TT TT OTT 

This can also be done in SAGE : 

SAGE 



sage: fl = lambda x: -1 
sage: f2 = lambda x: 2 

sage: f = Piecewise ( [ [ (0, pi/2 ) , f 1 ] , [ (pi/2 , pi ) , f 2 ] ] ) 
sage: PI = f.plotO 

sage: blO = [f . sine_series_coef ficient (n,pi) for n in range (1,10] 

sage: blO 

[2/pi, -6/pi, 2/(3*pi), 0, 2/(5*pi), -2/pi, 2/(7*pi), 0, 2/(9*pi: 
sage: sslO = sum ( [blO [n] *sin ( (n+1) *x) for n in range (len (b50) )] ) 

sage: sslO 

2*sin (9*x) / (9*pi) + 2*sin (7*x) / (7*pi) - 2*sin (6*x) /pi 
+ 2*sin (5*x) / (5*pi) + 2*sin (3*x) / (3*pi) - 6*sin (2*x) /pi + 2*sin(ic)/pi 
sage: b50 = [ f . sine_series_coef ficient (n, pi ) for n in range (1,50! 
sage: ss50 = sum ( [b50 [n] *sin ( (n+1) *x) for n in range (len (b) )] ) 
sage: P2 = sslO .plot (-5, 5, linestYle=" — ") 
sage: P3 = ss50 . plot (-5, 5, linestYle=" : " ) 
sage: (P1+P2+P3 ) . show () 



This illustrates how the series converges to the function. The 
function /(x), and some of the partial sums of its sine series, 
looks like Figure 4.6. 



Figure 4.6: f{x) and two sine series approximations. 



As you can see, taking more and more terms gives functions 
which better and better approximate f{x). 
The solution to the heat equation, therefore, is 

00 

u{x, t) = ^ bn{f) sin(— ) exp{-k{—)H). 

Next, we see how SAGE can plot the solution to the heat equa- 
tion (we use k = 1): 

I SAGE 

sage: t = var ("t") 

sage: solnSO = sum ( [b [n] *sin ( {n+1 ) *x) *e ^ ( - (n+1 ) " 2*t ) for n in range ( len (b50 ))] ) 

sage: solnSOa = sum { [b [n] * sin { (n + 1 ) *x) *e " ( - (n + 1 ) ^2* ( 1 / 1 ) ) for n in range ( len (b50 ))] ) 

sage: P4 = solnSOa . plot { , pi , linest yle=" : " ) 



sage: solnSOb = sum ( [b [n] *sin ( (n + 1) *x) *e" (- (n + 1 ) "2* ( 1/2 ) ) for n in range ( len (b50 ))] ) 

sage : P5 = soln50b.plot(0,pi) 

sage: solnSOc = sum ( [b [ n] *sin ( (n + 1 ) *x) *e " ( - (n + 1 ) " 2 * ( 1 / 1 ) ) for n in range ( len (b50 ))] ) 

sage : P6 = soln50c.plot(0,pi, linestyle=" -- " ) 

sage: (P1+P4 + P5+P6) .show() 



Taking 50 terms of this series, the graph of the solution at 
^ = 0, ^ = 0.5, ^ = 1, looks approximately like Figure 4.7. 




-0,5 



Figure 4.7: /(x), u{x,0.1), u{x,0.5), ^(x, 1.0) using 60 terms of the sine 
series. 



The heat equation with insulated ends boundary conditions 
models the temperature of an (insulated) wire of length L: 



d'^u{x,t) du{x,t) 

dx^ ~ dt 



Ux{0,t) = Ux{L,t) = 0. 



Here Ux{x^t) denotes the partial derivative of the temperature 
at a point x on the wire at time t. The initial temperature f{x) 
is specified by the equation 0) = f{x). 
Method: 

• Find the cosine series of f{x): 




n=l 



• The solution is 



oo 



U 



n=l 



Example: 

Let 




-1, 0<x<7r/2, 
2, 7r/2 < X < TT. 



Then L = tt and 




for n > and ag = 1. 
Thus 



/(^) y + cos(a:) + 02(7) cos(2a;) + 
This can also be done in SAGE : 

SAGE 



sage: fl = lambda x: -1 

sage: f2 = lambda x: 2 

sage: f = Piecewise ( [ [ (0, pi/2 ) , f 1 ] , [ (pi/2 , pi ) , f 2 ] ] ) 

sage: PI = f.plotO 

sage: alO = [ f . cosine_series_coef f icient (n, pi ) for n in range (10 

sage: alO 

[1, -6/pi, 0, 2/pi, 0, -6/(5*pi), 0, 6/(7*pi), 0, -2/(3*pi)] 

sage: a50 = [ f . cosine_series_coeff icient (n, pi ) for n in range (50 



sage: cslO = al0[0]/2 + sum ( [ al [n] *cos (n*x) 

sage: P2 = cslO.plot (-5, 5, linestYle=" — ") 

sage: cs50 = a50[0]/2 + sum ( [a50 [n] *cos (n*x) 

sage: P3 = cs50 .plot (-5, 5, linestyle=" : ") 

sage: (P1+P2+P3 ) . show ( ) 



for n in range (l,len 
for n in range (l,len 



] 

(alO) ) ] ) 
(a50) ) ] ) 



This illustrates how the series converges to the function. The 
piecewise constant function /(x), and some of the partial sums 
of its cosine series (one using 10 terms and one using 50 terms), 
looks like Figure 4.8. 

As you can see, taking more and more terms gives functions 
which better and better approximate f{x). 
The solution to the heat equation, therefore, is 

00 

u[x, ^) = y + 2^ ^nif) cos(-^) exp(-A;(— ) t). 

n=l 

Using SAGE , we can plot this function: 

I ^'^'^^ 1 

sage: solnSOa = a50[0]/2 + sum ( [a50 [n] *cos (n*x) *e " ( - (n+ 1 ) " 2 * { 1/100 ) ) for n in range ( 1 , len (a50 ))] ) 

sage: solnSOb = a50[0]/2 + sum( [a50 [n] *cos (n*x) *e" (- (n+1) "2* (1/10) ) for n in range ( 1 , len (a50 ))] ) 

sage: solnSOc = a50[0]/2 + sum( [aSO [n] *cos (n*x) *e" (- (n+1) "2* (1/2) ) for n in range (1, len (a50) )] ) 

sage: P4 = solnSOa.plot (0,pi) 

sage: P5 = soln50fa .plot (0, pi, linestyle=" : " ) 
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Figure 4.8: f{x) and two cosine series approximations. 



sage : P6 = soln50c.plot(0,pi, linestYle=" --" ) 
sage: (Pl+P 4+P5+P 6 ) .show() 



Taking only the first 50 terms of this series, the graph of the 
solution at t = 0, t = 0.01, t = 0.1,, t = 0.5, looks approximately 
like: 



-1 



2 



3 



Figure 4.9: f{x) = u{x,0), u(x,0.01), u{x,0.1), 0.5) using 50 terms of 
the cosine series. 

Explanation: 

Where does this solution come from? It comes from the method 
of separation of variables and the superposition principle. Here 
is a short explanation. We shall only discuss the "zero ends" 
case (the "insulated ends" case is similar). 

First, assume the solution to the PDE k^-^^ = ^"g^'^-* has the 
"factored" form 



u{x,t) = X{x)T{t), 
for some (unknown) functions X, T. If this function solves the 



PDE then it must satisfy kX"{x)T{t) = X{x)T'{t), or 



X"{x) _ 1 T'{t) 
X{x) ~ kl\tj' 

Since x, t are independent variables, these quotients must be 
constant. In other words, there must be a constant C such that 

Y(J)=kC, X"{x)-CX{x) = 0. 

Now we have reduced the problem of solving the one PDE to 
two ODEs (which is good), but with the price that we have 
introduced a constant which we don't know, namely C (which 
maybe isn't so good). The first ODE is easy to solve: 

T{t) = Aie^^\ 

for some constant Ai. To obtain physically meaningful solu- 
tions, we do not want the temperature of the wire to become 
unbounded as time increased (otherwise, the wire would simply 
melt eventually). Therefore, we may assume here that C < 0. 
It is best to analyse two cases now: 

Case C = 0: This implies X[x) = A2+Asx, for some constants 
A2,A^. Therefore 

u{x, t) = Ai{A2 + A^x) = Y + box, 

where (for reasons explained later) A1A2 has been renamed y 
and AiA^ has been renamed 6o- 

Case C < 0: Write (for convenience) C = — r^, for some r > 0. 
The ODE for X implies X{x) = A2C0s{rx) + ^3sin(ra;), for 
some constants A2,A3. Therefore 



u{x,t) = Aie '(^42 cos(rx)+743 sin(rx)) = (acos(rx)+6sin(rx))e 

where ^1^42 has been renamed a and ^1^43 has been renamed b. 

These are the solutions of the heat equation which can be writ- 
ten in factored form. By superposition, "the general solution" 
is a sum of these: 

t) = f^box^ ^^^li^n cos(r„x) + bn sin(rnx))e~^''"* 
= ^ + box + (ai cos(ri.T) + bi sin(ria:))e^^^i* 
+ (a2 cos(r2x) + 62 sin(r2x))e~^'^2* -)_ 

(4.5) 

for some a^, 5^, r^. We may order the r^'s to be strictly increasing 
if we like. 

We have not yet used the IC u{x, 0) = f{x) or the BCs u{0, t) = 
u{L, t) = 0. We do that next. 
What do the BCs tell us? Plugging in x = into (4.5) gives 

00 

= n(0, = ^ + E = T + + + - • 

z z 

n=l 

These exponential functions are linearly independent, so ao = 0, 
ai = 0, a2 = 0, ... . This implies 

u{x,t) = bpx+^y^^ bn sm{rnx)e~'^^^* = box+bism{rix)e''^'^^*+b2sm{r2x)e 

n=l 

Plugging in X = L into this gives 

= t) = boL + Y^ bn sin(r„L)e-^""^ 

n=l 



Again, exponential functions are linearly independent, so 60 = 0, 
bn sin (r^L) for n = 1, 2, .... In other to get a non- trivial solution 
to the PDE, we don't want bn = 0, so sin(rniv) = 0. This forces 
VnL to be a multiple of tt, say r„ = nir/L. This gives 



u{x,t) = y2^nsin{—x)e ''^"l)^* = i)^sm{—x))e ^^^^'^+62 sin(— a;))e 

^ 1j Lj Lj 

n=l 

(4.6) 

for some 6i's. The special case t = is the so-called "sine series" 
expansion of the initial temperature function u{x,0). This was 
discovered by Fourier. To solve the heat eqution, it remains to 
solve for the "sine series coefficients" hi. 

There is one remaining condition which our solution t) 
must satisfy. 

What does the IC tell us? Plugging t = into (4.6) gives 



f{x) = m(x, 0) = bn sin(— x) = 61 sin(— x))+62 sin(— x))+... . 

ij ij ij 

n=l 

In other words, if f{x) is given as a sum of these sine functions, 
or if we can somehow express f{x) as a sum of sine functions, 
then we can solve the heat equation. In fact there is a formula^ 
for these coefficients bn- 

bn = Y J /(^) cos{^x)dx. 
It is this formula which is used in the solutions above. 



Exercise: Solve the heat equation 

^Fourier did not know this formula at the time; it was discovered later by Dirichlet. 



2^d u{x,t) du{x,t) 

dx^ dt 



UxiO.t) = Ux{3,t) = 
u{x, 0) = X, 

using SAGE to plot approximations as above. 



4.4 The wave equation in one dimension 

The theory of the vibrating string touches on musical theory 
and the theory of oscillating waves, so has likely been a concern 
of scholars since ancient times. Nevertheless, it wasn't until 
the late 1700s that mathematical progress was made. Though 
the problem of describing mathematically a vibrating string re- 
quires no calculus, the solution does. With the advent of cal- 
culus, Jean le Rond dAlembert, Daniel Bernoulli, Leonard Eu- 
ler, Joseph-Louis Lagrange were able to arrive at solutions to 
the one-dimensional wave equation in the eighteenth-century. 
Daniel Bernoulli's solution dealt with an infinite series of sines 
and cosines (derived from what we now call a "Fourier series", 
though it predates it), his contemporaries did not believe that 
he was correct. BernouUis technique would be later used by 
Joseph Fourier when he solved the thermodynamic heat equa- 
tion in 1807. It is Bernoulli's idea which we discuss here as well. 
Euler was wrong: Bernoulli's method was basically correct after 
all. 

Now, d'Alembert was mentioned in the lecture on the trans- 
port equation and it is worthwhile very briefly discussing what 
his basic idea was. The theorem of dAlembert on the solution 
to the wave equation is stated roughly as follows: The partial 
differential equation: 

= c • 

is satisfied by any function of the form w = w{x,t) = g{x + 
ct) + h{x — ct), where g and h are "arbitrary" functions. (This is 
called "the dAlembert solution".) Geometrically speaking, the 



idea of the proof is to observe that ^ ± is a constant times 
the directional derivative Dy-^w{x,t), where v± is a unit vector 
in the direction (±c, 1). Therefore, you integrate 

Dy-Dy-^w{x,t) = (const. - c • = 

twice, once in the direction, once in the vl, to get the solu- 
tion. Easier said than done, but still, that's the idea. 

The wave equation with zero ends boundary conditions models 
the motion of a (perfectly elastic) guitar string of length L: 

2 d'^w{x,t) _ d'^w{x,t) 
dx^ dip" 

w{0,t) = w{L,t) = 0. 

Here w{x,t) denotes the displacement from rest of a point x on 
the string at time t. The initial displacement f{x) and initial 
velocity g{x) at specified by the equations 

w{x,0) = f{x), Wtix.O) = g{x). 

Method: 

• Find the sine series of /(x) and g{x): 

oo 00 

fix) ~ X]&n(/)sin(— ), g{x) ~ J^^M sm{—). 

n=l n=l 

• The solution is 

/ X / r\ / nTTt. LbJg) . , nnt.. . .mix. 

w{x,t) = > (6„(/)cos(c— )+ sm(c— ))sm(— — ). 

L cmr L L 

n=l 



Example: Let 



1, 0<t<7r/2, 



2, 7r/2 <t < TT, 
and let g{x) = 0. Then L = tt, bn{g) = 0, and 

, /.N 2 r , . , ^ 2 cosfnTr) - 3 cos(l/2n7r) + 1 
= - / f{x) sm{nx)dx = -2 ^ — ^— ^ ' ' 



TT 



n 



Thus 



2 6 2 
/(x) ~ &i(/) sin(a:)+62(/) sin(2a:)+... = — sin(a:) sin(2a:)+— sin(3a:)+. 

The function /(x), and some of the partial sums of its sine series, 
looks like 



-3 -2 



Figure 4.10: Using 50 terms of the sine series of f{x). 



This was computed using the following SAGE commands: 

SAGE 



sage: fl = lambda x: -1 
sage: f2 = lambda x: 2 



sage: f = Piecewise ( [ [ ( , pi/2 ) , f 1 ] , [ (pi/2 , pi ) , f 2 ] ] ) 

sage: PI = f .plot (rgbcolor= (1, 0, 0) ) 

sage: b50 = [ f . sine_series_coef f icient (n, pi ) for n in range (1,50! ] 

sage: ss50 = sum ( [b50 [i-1 ] *sin (i*x) for 1 in range (1, 50) ] ) 

sage: b50[0:5] 

[2/pi, -6/pi, 2/(3*pi), 0, 2/(5*pi)] 

sage: P2 = ss50 .plot (-5, 5, linestyle=" — ") 

sage: (P1+P2) . show ( ) 



As you can see, taking more and more terms gives functions 
which better and better approximate f{x). 
The solution to the wave equation, therefore, is 

w{x, t) = y{bn{f) cos(c— ) + sm(c— )) sm(— — ). 

^ — ^ L cniT L L 

n=l 

Taking only the first 50 terms of this series, the graph of the 
solution at t = 0, t = 0.1, t = 1/5, t = 1/4, looks approximately 
like: 



^ ■ 



Figure 4.11: Wave equation with c = 3. 
This was produced using the SAGE commands: 



SAGE 



sage: t = var ("t") 

sage: w50tl = sum ( [b50 [ i-1 ] * sin (i*x) *cos ( 3 *i * { 1/ 10 ) ) for i in range (1, 50) ] ) 

sage : P3 = w5 Ot 1 . plot ( , pi , linestyle=" : " ) 

sage: w50t2 = sum ( [b50 [ i-1 ] * sin (i*x) *cos ( 3 *i * ( 1/5 ) ) for i in range ( 1 , 50 )] ) 

sage: P4 = w50t 2 . plot ( , pi , linestyle=" : " , rgbcolor= { , 1 , ) ) 

sage : w50t 3 = sum ( [b50 [i-1] *sin(i*x) *cos (3*i* (1/4) ) for i in range (1,50) ] ) 

sage: P5 = w50t 3 . plot ( , pi , linest yle=" : " , rgbcolor= ( 1/3 , 1 / 3 , 1 /3 ) ) 

sage: (P1+P2+P3+P4+P5 ) .show() 



Of course, taking terms would give a better approximation to 
w{x^t). Taking the first 100 terms of this series (but with dif- 
ferent times): 




Figure 4.12: Wave equation with c = 3. 



Exercise: Solve the wave equation 



i^d'^wix.t) d'^w{x,t) 

w{Q,t) = w{3,t) = 
w{x, 0) = X 

using SAGE to plot approximations as above. 
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